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Abstract 

We present a numerical model for the dynamics of thin viscous threads based 
on a discrete, Lagraugian formulation of the smooth equations. The model 
makes use of a condensed set of coordinates, called the centerline/spin repre- 
sentation: the kinematical constraints linking the centerline's tangent to the 
orientation of the material frame is used to eliminate two out of three degrees 
of freedom associated with rotations. Based on a description of twist inspired 
from discrete differential geometry and from variational principles, we build a 
full-fledged discrete viscous thread model, which includes in particular a dis- 
crete representation of the internal viscous stress. Consistency of the discrete 
model with the classical, smooth equations is established formally in the limit 
of a vanishing discretization length. The discrete models lends itself naturally 
to numerical implementation. Our numerical method is validated against refer- 
ence solutions for steady coiling. The method makes it possible to simulate the 
unsteady behavior of thin viscous jets in a robust and efficient way, including 
the combined effects of inertia, stretching, bending, twisting, large rotations and 
surface tension. 

Keywords: viscous rod, bending, twisting, Rayleigh- Taylor analogy, viscous 
coiling 
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1. Introduction 

1.1. Context 

The flow of thin viscous filaments is relevant to a variety of industrial pro- 
cesses such as the drawing and spinning of polymer and glass fibers [3 HI 13] , and 
to natural phenomena such as formation of Pele's hair by lava ejected at high 
speed by volcanoes [1]. In art, Jackson Pollock took advantage of the coiling 
instability of a thin viscous fluid, the paint, impinging a surface, the canvas, 
to produce a variety of decorative patterns by a fluid- mechanical process [5]. 
A commonplace version of the same coiling instability is observed when a thin 
thread of honey is poured on a morning's toast. This steady coiling problem is 
prototypical of the dynamics of thin threads. Its apparent simplicity has made 
it appealing to fluid mechanicians for a long time [6l [7] ; however the various 
regimes of steady coiling and its non-linear features, such as multistability, have 
been understood in full details only recently [5J [51 UHl E] • To a large extent, the 
analysis of steady coiling has been made possible by the availability of numerical 
simulation: the shape of the thread in the co-rotating frame is stationary, and is 
given by a non-linear boundary-value problem [5] which has been solved using 
numerical continuation. 

In this paper we are interested in the simulation of the unsteady behavior of 
thin threads, which is far less advanced. As an illustration, consider a recently 
proposed variant of the coiling problem, similar to Pollock's painting technique, 
whereby the target surface moves horizontally at a constant velocity. The rel- 
ative motion suppresses steady coiling solutions and forces the flow to become 
unsteady. More than ten different patterns can be produced by varying the 
lateral velocity of the surface and the fall height [HI [T3j , a number of which 
have convoluted and intriguing shapes. The patterns are reminiscent of stitch 
patterns, and the experiment has been coined the 'fluid-mechanical sewing ma- 
chine'. This experiment nicely illustrates the complex behavior that can result 
from the dynamics of a thin, perfectly viscous fllament. Existing numerical 
methods are unable to reproduce this behavior, even though the principle of the 
experiment is simple. This highlights the need for a robust and efficient method 
for simulating the dynamics of thin viscous threads. 

The dynamics of thin viscous threads is governed by the interplay of three 
local modes of deformation, namely stretching, bending and twisting modes [141 
(TSj. At the global scale, these modes are coupled by geometrically non-linear 
terms accounting for finite rotations. This coupling makes the resulting dy- 
namics remarkably rich. Another, unfortunate consequence of the nonlinearity 
is the absence of analytical solutions to the dynamical equations. This makes 
the development of robust simulation methods even more desirable. The main 
difficulty in developing such a method is that the underlying equations are 
numerically stiflF, as the governing equations are non-linear partial differential 
equations of fourth order in space. This paper tackles this difficulty by intro- 
ducing a careful and well-controlled space discretization. In fact, we introduce a 
full-fledged discrete viscous thread model by extending all the relevant physical 
quantities, such as strain rates and internal stress, to the discrete setting. 
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Fluid mechanical problems involving free boundary conditions can be simu- 
lated using refined variants of the marker and cells method, namely the method 
of O Hi] for 2d viscous flows, and the GENSMAC method [El [IS] for 3d 
viscous flows; more recently implicit schemes coupled with projection methods 
have been proposed, see e. g. HHj ■ The present paper is concerned with thin fila- 
ments, for which the above methods are not efficient: when the thickness is small 
compared to the longitudinal length scale, it is beneficial to use dimensionally 
reduced equations as a starting point for simulations. Thanks to dimensional 
reduction, the structure of the flow at small scale is solved analytically, which 
makes it possible to use a simulation grid much coarser than the thickness. 

While our simulation method addresses the general non-steady dynamics of 
a thin thread governed by the combined effects of twist, bending and stretch- 
ing forces, inertia and large rotations, a number of particular cases have been 
simulated in the literature. Steadily rotating viscous threads are described by 
time-independent equations in the co-rotating frame, which have been solved 
numerically using methods for two-point boundary value problems [51 llli 121] . 
The dynamics of a viscous string, where both the bending and twisting modes 
are neglected, has been considered [52] . The periodic folding of a viscous thread 
or sheet has been considered in a 2d geometry j23l El] where twist does not 
play any role. By combining the simulation of steady solutions with analytical 
expansions describing oscillatory perturbations of small amplitude, the stabil- 
ity of both the steady coiling solution [TT] and of the catenary-like profile of a 
dragged thread have been calculated. Many other problems, such as the 
existence of rotatory folding, the competition between folding and coiling [26] . 
the stitch patterns produced by the fluid-mechanical sewing machine [12j and 
the destabilization of steady coiling by precession [27] remain inaccessible to 
those simulation methods that are based on restrictive assumptions. 

In comparison to viscous threads, elastic rods have received a lot of attention, 
both from the perspective of analysis [281 ISS [SQl E] and simulation [32l [33l 
[34l [35l [36l l37] . By the Rayleigh- Taylor analogy [7], the stress in a viscous 
fluid is identical the stress in an elastic solid having the same geometry, when 
the strain rate relevant to the viscous case is replaced with the current strain 
relevant to the elastic case. Stated differently, the main difference between the 
elastic and viscous problems is the presence of an additional time derivative 
in the right-hand side of the viscous constitutive laws. This analogy explains 
the buckling of viscous sheets [351 133] : a phenomenon usually associated with 
elastic structures. One can take advantage of this analogy to simulate the 
dynamics of viscous threads using a simulation tool written for elastic rods |40] ; 
we explored this approach in a conference paper [ ¥T] . The possibility to recycle 
an existing elastic code for viscous simulations with minimal additional work is 
very attractive. However, the initial effort associated with implementing and 
validating an elastic code is high. In situations where no elastic code is available, 
it is simpler to implement a viscous simulator directly. In the present paper we 
follow this approach and propose a numerical method that does not rely on an 
external library for solving the dynamics of elastic rods. 

While thin elastic rods are usually considered inextensible, the stretching 



3 



mode has to be retained in viscous threads, in addition to the usual bending 
and twisting modes. In both the elastic and viscous cases, dimensional analysis 
shows that the strain, or strain rate, associated with the stretching mode is 
small compared to that associated with the bending and twisting modes, as 
the corresponding modulus is larger by a factor proportional to the inverse of 
the small aspect-ratio squared, a very large number. A specificity of the viscous 
case is that the stretching mode cannot be neglected, even though its strain rate 
is small. This is well illustrated by the phenomenon of helical coiling: a thin 
thread poured from a container onto a fixed obstacle gets stretched by gravity 
and remains straight over almost the entire fall height, but it bends and twist 
severely in a small boundary layer near the bottom. Even though the stretching 
is very mild, its effect is cumulated over the entire time of descent, unlike the 
bending and twisting modes that come into play only near the surface. For 
this reason, the stretching mode has to be considered in simulations of viscous 
threads. By the incompressibility condition, variations of the thread's radius 
along its centerline need be considered as well. Another difference with the 
elastic case is that capillary forces can have a strong effect on the motion, and 
must be taken into account. 

1.2. Model 

The derivation of dimensionally reduced models for thin viscous filaments has 
a long history and is still a research topic. The equations for thin viscous threads 
were derived by asymptotic expansion from the equations for a 3d viscous flu- 
ids by Entov and Yarin [42j. Their work builds upon the previous analyses of 
viscous stretching by Trouton [T3], and of viscous bending by Buckmaster and 
co-workers [43l|44]. In the case of elliptical cross-sections, the dynamics of the 
centerline and the evolution of the geometry of the cross-section are coupled; the 
corresponding equations have been derived by dimensional reduction in Id by 
Dewynne and collaborators [35^, and later extended to a capillary fluid by |3S]. 
Recent derivations of the equations for thin threads beneflt from a clear iden- 
tification of the mechanical quantities in the Id model model [17], and of the 
systematic use of Lagrangian coordinates |21j . Asymptotic models accounting 
for more general constitutive laws have been proposed: the case of a visco-elastic 
fluid is treated in [48j , and a general framework is considered in which can 
produce a variety of asymptotic models when a specific set of physical effects is 
considered. 

Here we consider the dynamics of a thin filament of an incompressible, purely 
viscous fiuid having circular cross-section, under the action of external forces 
such as gravity, and internal forces (viscous stretching, bending, twisting, and 
capillary tension). We consider the 3d problem, and the curvature and kine- 
matical twist of the thread can be comparable to, or smaller than the inverse of 
the thread's length. Even though the fluid is very viscous, the effect of inertia is 
considered. The role of inertia is well illustrated by the classical analysis of the 
pendulum modes of a viscous string, see e. g. [TUl: in this almost straight ge- 
ometry, the flow in the axial direction is typically governed by a small Reynolds 
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number and dominated by viscosity, although the flow in the transverse direc- 
tion, which is characterized by different length and time scales, is associated 
with a much larger Reynolds number and dominated by inertia. In the general 
case, the local axial and transverse modes get coupled by the curvature of the fil- 
ament. For this reason, we retain the inertial term in the balance of momentum 
even though the constitutive laws are dominated by viscosity (Stokes' fluid). 

We assume that the cross-section of the filament is a disk, and remains 
so as the filament deforms. Even when a viscous filament is extruded from a 
non-circular opening, surface tension tends to round off the shape of the cross- 
sections; this happens over a time scale which we assume is short compared 
that of the flow. This is not always the case, and the possibility to account for 
non-axisymmetric cross-sections has been demonstrated in the elastic case using 
our numerical method [35l [41] . While the general case raises no fundamental 
difficulty, our presentation is limited to the case of axisymmetric cross-sections, 
which is simpler: all the directions in the cross-sections are then equivalent and 
there is no need to keep track of their absolute orientations. The case of tubes, 
i. e. of non simply-connected cross-sections [S^, is not considered here but this 
extension may be considered in future work as well. 

1.3. Proposed approach 

The main features of our numerical method are the following. It is based on 
a Id model obtained by dimensional reduction, which makes it much more effi- 
cient than a general-purpose model for 3d viscous flows. We use a Lagrangian 
grid, making simulation vertices flow along with the fluid; this simplifies the 
computation of viscous forces which, by the constitutive law, are proportional 
to the comoving time derivative of the kinematical twist and curvature. We 
use a reduced set of coordinates and eliminate two out of three degrees of free- 
dom associated with rotations, by making use of the fact that cross-sections 
initially perpendicular to the centerline remain so during motion: rotations are 
represented using a single degree of freedom. In addition, the absolute angle of 
twist of the cross-sections is eliminated from the equations, using the fact that 
the cross-sections are circular. This results in an effective description of rota- 
tions which only makes use of the instantaneous angular twist velocity, denoted 
V. Internal viscous forces include the three physically relevant contributions 
of stretching, bending and twisting. The discrete expression of these forces is 
derived based on a Rayleigh potential using variational methods. This leads 
naturally to discrete viscous forces. The Rayleigh potential plays a similar role 
in the viscous setting as the elastic energy in the elastic setting. The viscous 
forces are computed in a linear implicit manner, which means that the expres- 
sion for the viscous forces is extrapolated to their value at the end of the time 
step, based on a linearization carried out at the beginning of the time step. 
This provides a good compromise between stability and ease of implementa- 
tion. A fully implicit evaluation of these forces is also possible, as discussed in 
section 19.31 

Our main contribution is to propose an elaborate space discretization of a 
viscous thread by considering a fully discrete model. All the relevant physical 
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quantities, such as the strain rates and the viscous forces, are defined in the 
discrete setting. In this discrete view, bending is represented by the turning 
angles of the centerhne at the vertices, and we do not need to assume that the 
turning angles remain small at all times. As a consequence, the simulation of 
the discrete model is stable even if the mesh size is comparable to the smallest 
radius of curvature of the thread. By contrast, stability of a numerical scheme 
based on an ad hoc discretization of the smooth equations usual requires that 
the grid size is much smaller than the smallest radius of curvature; in practice, 
this requirement is severe as thin threads tend to spontaneously form strongly 
curved region near the endpoints (such as the tiny rotating coils at the bottom 
of a thin thread stretched by gravity and hitting a hard surface). The ability 
to run simulations with a relatively coarse mesh is an important advantage of 
the discrete model, as the maximum time step allowed by Nyquist stability 
criterion decreases rapidly with mesh size in the presence of fourth-order space 
derivatives: the simulation of coarser meshes is much more efficient. Given the 
numerical stiffness of the underlying problem, robustness is a central issue in 
the simulation of viscous threads. We address this issue by keeping full control 
of the space discretization. 

The discretization of bending in elastic rods is routinely done using flexural 
springs at hinges [51], and the extension to viscous bending is straightforward. 
Enforcing the twist forces is much less common. Our discretization of twist is 
based on a discrete notion of twist, which is directly borrowed from our previous 
work on elastic rods [3S] . This discrete notion of twist is based on concepts from 
discrete differential geometry, namely the holonomy of a discrete curve. This 
representation of twist builds upon previous work highlighting the geometrical 
origin of twist [521 1531 154] , and on related numerical methods used in mechanical 
engineering |55| . 

This paper is organized as follows. In section [2] we derive useful identities of 
geometry and differential calculus. In section [3] the equations for thin viscous 
threads are formulated in a way that prepares the extension to the discrete the 
setting, by making use of a Lagrangian description of motion and by deriving 
the internal viscous forces and moments from variational principles. Section [4] 
establishes the equivalence with the formulation of the equations for thin threads 
classically used by fluid mechanicians. In section [5] the discrete model is pre- 
sented in close analogy with the smooth case of section [3] this section is at the 
core of the paper. Section [6] considers time discretization; formulae that are 
required in the implementation are summarized, and we discuss the treatment 
of boundary conditions and adaptive mesh refinement. In section [7| we consider 
the coupling of the thread with external bodies, such as a fluid container or 
hard surface. In section |8l the code is validated against reference solutions for 
steady coiling. In section]9j we discuss limitations and perspectives. 

In this paper, an effort is made to establish the equivalence of different 
formulations. The formulation of the equation of motion due to Kirchhoff is 
popular among fluid mechanicians as it has a clear intuitive meaning, while 
that based on the Rayleigh potential is rarely if ever used in this community 
but lends itself to a natural discretization. Working out the connection between 
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these equivalent formulations will hopefully help to make this paper accessible 
to different communities. In addition, presenting the discrete model in parallel 
with the smooth model provides mutual insights into them. The reader should 
keep in mind that only a small fraction of the formulae presented in the paper 
are required for the purpose of implementing the method. These formulae are 
recapitulated in section |6.4[ 



2. Mathematical toolbox 



In this section we introduce some notations and mathematical identities 
relevant to the mechanics of thin viscous threads, such as infinitesimal rotations 
and the covariant derivative. 



2.1. Perpendicular projection 

We use underlines for vectors, and double underlines for rank- two tensors 
(matrices). For any unit vector q and for any vector a, the projection Pj^ in the 
direction perpendicular to q is defined by 

P±{q,a) ^ {I- q(E) q) ■ a- {q- a)q, (1) 

where the last term in the right-hand side is the longitudinal projection. 



2.2. Infinitesimal rotations 

Consider an orthonormal frame d^{a), i ~ 1,2,3, which is function of a 
continuous parameter a: for any pair of integer indices 1 < «,J < 3, and for 
every cr, we have 

dd'y) ■ d^{<j) = (2) 

where Sij is Kronecker's symbol, equal to 1 ii i — j and to otherwise. We shall 
assume that the frame is smooth with respect to the parameter a. 

Infinitesimal rotations can be expressed by means of a skew-symmetric ma- 
trix, or cquivalently as the multiplication by a vector T{a) using the cross prod- 
uct. More accurately, for any value of the parameter cr, the Darboux vector 
r(t7) is defined as the unique vector such that for any i — 1,2, 3: 

^=r(a)xd». (3) 

This definition will be used both when cr is the time t, or the arc length S. 
This leads to the notions of instantaneous angular velocity £(<) = lo, or twist- 
curvature vector r(5) — tt, respectively. 

An explicit expression for the Darboux vector can be found by singling out 
any particular vector in the triad, say d^: 

m = d,x^ + r,d„ (4) 
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where T3 = F • dg is defined by 

r3(a) = ^.d,(.). (5) 

This can be checked by inserting equation ^ into equation Q and then into 
equation ([s]). 

2.3. Covariant derivative 

Let us consider an orthonormal triad ^^(cr) and the associated Darboux 
vector r(cr). For any vector field a(cr), we define the covariant derivative as 

%M^%M_r(a)xa(a). (6) 
da do- 
lt can be interpreted as the derivative measured in the frame moving with the 
triad d^. In agreement with this interpretation, we note that equation ([s]) can 
be rewritten as 

= 0. (7) 

da 

The covariant derivative satisfies the following identities, the proof of which 
is left to the reader. The general Leibniz rule for the product of a scalar function 
/(cr) by a vector a(cr) reads 

m^M^ = a{a) + f{a)~^. (8) 

dcr dcr do- 

In the first term of the right-hand side above, the regular derivative of a scalar 
function is used. In the case of the scalar product of two vectors we have 

d(a ■ b) da , db 

, -' = -b + a- ^. (9) 
do- da da 

In the context of thin viscous threads, an important property of the covariant 
derivative is its compatibility with the tangent and the normal projections. For 
any vector a we define the tangent and normal projections, using the last triad 
vector dg as the tangent direction: 

a = a3C?3 + a^, where ag = a • dg and a^ = (dg , a) . (10) 

Here Pj^ is the perpendicular projection operator defined in equation ([T]). A 
consequence of equations ([7]-[9]) is that the covariant derivative lives in the plane 
orthogonal to dg if a{a) does so for all values of a, and is aligned with the 
tangent if a(cr) is aligned with the tangent for all a. In other words, both the 
tangential and normal projections commute with the covariant time derivative, 

~^.d,=^, (~^Y=~^. (11) 

da da yda J da ^ ^ 

By contrast the regular derivative d/da does not commute with the projection 
operators. Note that in the first equation above, we use a regular derivative for 
the scalar quantity ag, hence the absence of a tilde. 
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3. Smooth setting: a Lagrangian description of viscous threads 



The equations for the dynamics of thin viscous threads are usually expressed 
in Eulerian variables, see for instance references [111. In the present section we 
reformulate these equations in Lagrangian variables. This will make it possible 
to extend the geometrical discretization of twist of Bergou et al. |35j , and then 
to derive a discrete model of viscous threads in a natural way. A Lagrangian 
formulation of the equations can be found in reference j21| but we are not aware 
of any numerical method that actually makes use of it; most numerical papers 
consider steady problems, for which an Eulerian grid is sufhcient. 

Below, we introduce the smooth quantities that are relevant to the kine- 
matics and dynamics of viscous threads. They are introduced in a way that 
makes it straightforward to identify their discrete counterparts later on. Since 
the viscous thread can stretch, we make a careful distinction between the arc 
length measured in reference configuration, which is denoted S, and the arc 
length measured in actual configuration, denoted s. 

3.1. Reference configuration 

Our Lagrangian description requires the definition of a reference configura- 
tion. This reference configuration can be imaginary, and does not need to be 
the configuration of the thread at any particular time. A convenient choice is to 
define the reference configuration to be an infinite, circular cylinder of constant 
radius oq, as illustrated in figure [l] Obviously the equations of motion will be 
independent of the choice of the reference configuration, and of the value of the 
reference radius oq in particular. 

The fluid is considered incompressible. It therefore simpler, although it is 
not required strictly speaking, to assume that the mapping between the refer- 
ence and actual configurations preserves volume. This is what we do in equa- 



tions (23 1. 



3.8. Kinematics of centerline 

Let S be the Lagrangian coordinate, and t the time. For any function f{S, t), 
we denote its spatial derivative using a prime, 

m<) = «^. (12) 

and its time derivative using a dot, 

M0 = ^. (13) 

By definition, the Lagrangian coordinate S follows fixed fluid particles. The 
time derivative introduced above is known as a convective derivative, and is 
often written f — ^ — the Eulerian context. 

At time t the centerline of the thread is given by the function x{S,t), see 
figure [ij The material tangent of the thread is denoted T{S,t) and defined by 
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actual 

Figure 1: Reference and actual configuration of the thread. 



T{S,t)=x'{S,t). 



(14) 



Note that this is not a unit vector in general. 

The norm of T^{S,t), denoted £{S,t), measures the amount of stretching of 
the centerhne with respect to the reference configuration: 



£iS,t) = \T{S,t)\ 



dx{S, t) 



dS 



The unit tangent to the centerhne is then defined as 

T{S,t) 



t{S,t) 



£{S,t) 



(15) 



(16) 



In our Lagrangian description, the arc length s in actual configuration is 
viewed as a secondary quantity. It can be reconstructed by integration of the 
differential equation expressing the identity ds — \dx\, namely 



The Lagrangian axial strain rate d is defined by 

d£{S, t) 



d{S,t) 



dt 



(17) 



and measures the rate of stretching per unit time. Note that this quanti ty di ffers 
from the Eulerian strain rate, noted and defined later in equation (65c I. 



In the Lagrangian framework the velocity u is simply the time derivative of 
position, 

dx{S, t) 



u{S,t) 



dt 



(18) 
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and the acceleration its second time derivative, 

.... du{S,t) d\{S,t) 

A kinematical relation between the strain rate d{S, t) and the velocity u{S, t) 
can be established as follows: 

"'^^'^^^ at ~ 2^ dt - 2i dt -I- dt--^^'*> ds • ^^^> 

In the last equality, we have used the following identity, coming from the per- 
mutation of derivatives with respect of t and S: 



dT 


-J 


(dx\ 




fdx\ 


du 


dt 


dt 


KdS 


dS 


Kdt 


dS 



Another useful identity follows from inserting the definition of £ in equa- 



tion (16), T — it, into the left-hand side of the equation above. Expanding the 
time derivative, we find £t + £i = u'{S). Applying the perpendicular projection 
operator Pj_{t, •) on both sides and using Pj^{t,t) = and Pj_{t,i) = i (the 
other term cancels since t-t=^ ^ = 0) i we have 



dt £{S,t) y-'^ ^ " dS 

This equation yields the time derivative of unit tangent as a function of center- 
line geometry and of the velocity u. 

3.3. Incompressibility: radius and related quantities 

The viscous fluid is considered incompressible, as explained in section |3.1[ 
and we assume that the mapping from the reference configuration to the actual 
configuration preserves volume. The reference configuration is a cylinder with 
a uniform radius oq and a cross sectional area Aq = ttuq. The radius of the 
thread in the actual configuration is denoted a{S,t), the area is A{S,t), and 
I{S, t) stands for the moment of inertia. We assume that the thread has locally 
a cylindrical geometry, in which case 

A{S,t)=7Ta\S,t), J(5,i) = ^^^^. (22) 

The fluid volume enclosed in an infinitesimal length of the thread is ^(5", t) ds = 
A{S, t) £{S, t) dS in the actual configuration, and Ao{S, t) dS in reference config- 
uration, see figure [ij As a result, the incompressibility of the fiuid is expressed 

by 

where the subscript naught refers to the reference configuration, for which we 
have ^0 = 1 by convention. In particular, the moment of inertia in the reference 
configuration reads /q = ira'^/A:. 
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I dt 




Figure 2: Darboux vectors (infinitesimal rotations) associated with spatial or temporal deriva- 
tive of the material frame: instantaneous rotation veloci ty a; and twist-cur vatu re vector tt. A 
compatibility condition for and tt is derived in section[3.5[ see equation (jSSl. 



3.4^. Material frame, angular velocity 

To complete the description of the motion, we need to keep track of twist, 
defined as the rotation of the cross sections about the tangent. Indeed twist 
gives rise to viscous shear stress in the plane of a cross-section, which affects 
the dynamics of the thread. With the aim to measure twist, we introduce a 
triad, denoted {di{S,t),d2{S,t),(U,{S,t)), which is rigidly attached to the cross- 
sections. This triad follows the motion of the surrounding fluid particles and is 
called the material frame. 

The flow inside the thread is shearless in the limit of very thin thread, as 
explained for instance in the work of Ribe and collaborators [TT]. A similar 
argument holds in the case of elastic rods, see for instance reference [S5]. As a 
result, material cross sections cannot slide with respect to another and remain 
perpendicular to the local tangent to the centerline. This is known as Kirchhoff 
kinematical hypothesis, but it can be justifled rigorously by asymptotic anal- 
ysis. Therefore we impose the following conditions: (?) the triad satisfies the 
orthonormality condition ([2]), and [ii) it is compatible with the centerline in the 
sense that 

d^{S,t)=t_{S,t). (24) 

This kinematical condition is at the core of the mechanics of thin threads; it 
couples the rotations of the material frame, measured using the triad dj in the 
left-hand side of the equation, with the motion of the centerline whose tangent 
appears in the right-hand side. 

The definition of the Darboux vector in equation ([s]) yields, after replacing 
the general parameter a with time, a = t, 

^^^^u,{S,t)xd,{S,t) (25) 

for any value of the index i — 1,2,3. In the context of time derivation, the 
Darboux vector T{t) is denoted oj{S, t) and is called the angular velocity vector; 
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its tangential component is called the axial spin velocity Fa — v{S,t), 



v{S,t) ^ui{S,t) ■t{S,t). (26) 

By equation ([H]), u is given by f = di'd2- explicit expression for the angular 
velocity is obtained from equation Q: 

ui{S, t) = t_{S, t) X t_{S, t) + v{S, t) t{S, t). (27) 

A second Darboux vector is obtained in the case of spatial derivatives, insert- 
ing a — S iw equation This vector is associated with infinitesimal changes 
of the Lagrangian coordinate S, and is called the twist-curvature vector. It is 
denoted 7r{S,t) and satisfies the fundamental relation 

^^^^^=7L{S,t)xd,iS,t), (28) 

for any value of the index 1 < « < 3. Its tangential component is called the 
(Lagrangian) kinematical twist and noted t{S, t) instead of the generic notation 

T{S,t)^TT_{S,t)-t{S,t). (29) 

This scalar T{S,t) measures the rate of rotation of the material frame about 
the tangent: T{s,t) = d!i ■ d^- It is different from the Frenet-Serret notion of 
torsion for a three-dimensional curve which does not make any reference to 
the material frame d^. The Frenet-Serret torsion measures the non-planarity 
of a curve, although the kinematical twist r can be non-zero even though the 
centerline is planar, as happens for instance in the case of a straight but twisted 
configuration. 

The normal projection of tt is given by the general expression of the Darboux 
vector in equation Q, combined with the condition of compatibility (24), 



TL{S,t)^ K{S,t) + T{S,t)t_{S,t), (30) 
where we have introduced the Lagrangian binormal curvature 

^(5,t)=t(5,0x (31) 

Note K{S, t) depends only on the centerline, not on the material frame. Con- 
sistently with our Lagrangian approach, both the kinematical twist r and the 
binormal curvature refer to a unit increment of the arc length S in the ref- 
erence configuration, not to the arc length s in the actual configuration. As a 
result, these quantities differ from the usual twist and curvature used in the Eu- 
lerian framework and some care is required when comparing our equations with 



their Eulerian variants. Equation (30) shows that the tangential component of 
TT encodes the twist while its normal projection encodes for the curvature, hence 
the name twist-curvature vector. 
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The measure of strain relevant to the stretching mode is provided by the 



extension £{S,t) defined earlier in equation (16). This vector Tr{S,t) just intro- 
duced provides a measure of strain relevant to the two other modes, twisting 
and bending. For viscous threads, we need to evaluate the rates of strain. This 
is the goal of the next sections: in section |3.5| we derive an identity for the 



covariant time derivative of Tr{S,t), and in section 3.6 we identify this quantity 



as the rate of strain, which is needed in the constitutive relations. 



3. 5. Compatibility of rotations 

The rotation velocity u and the twist-curvature vector tt characterize in- 
finitesimal changes of the material frame d^{S,t) corresponding to increments 
of time t and of arc length S, respectively. A compatibility condition relating 
the vectors u/ and can be derived from the identity of the cross-derivatives in 
Lagrangian coordinates {t,S): 



d_ f dd,it,s) 

dt V dS 



d_ f dd,{t, S) 
dS [ dt 



In this expression, we express the innermost derivatives with the help of equa- 



tions ( 25 1 and ( 28 ) and write 

d{n X dj) d{uj x d^) 



dt 



dS 



We use equations (25) and (28 1 again to expand the derivatives, and find 



dn 
'dt 



dS 



X = — TT X (w X dj) + cj X (tt X d^). 



The right-hand side can be simplified using Jacobi's identity, ax {bxc) + bx 
(c X a) + c X (a X 6) = 0, valid for any set of vectors a, b, c. With a = tt, 6 = cj 
and c = d^, this yields 

which is also known as the Maurer-Cartan identity. Here we have simplified by 
dj, since the relations hold for any value of the index i = 1, 2, 3. 

In the right-hand side we can identify the covariant derivative defined in 
section [2.3[ here with a = t and t — d^- It is defined for an arbitrary vector 
field a{S, t) by 

da{S, t) _ da{S, t) 
dt ~ dt 

We can then rewrite the compatibility condition in a compact form, 



-ui{S,t) X a{S,t). 



(32) 



doj{S, t) dniS, t) 



dS 



dt 



(33) 
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Note that the role of lo and tt is symmetric and it is possible to rewrite this 
equation the other way around, 

d7L{S, t) _ doj{S, t) 



dS 



dt 



However, we shall only make use of equation (33) in the following, as it provides 
the relevant measure of rate of strain. 



3.6. Strain rate vector 



The quantities appearing in equation ( 33 ) are fundamental for the dynamics 
of viscous rods, and are related to the strain rates of the different modes of 
deformation. The left-hand side is denoted e: 



dui{S, t) 
dS ' 



(34) 



The vector e, called the strain rate vector^ provides the measure of the rates of 
deformation relevant to the twist and curvature modes. As such it appears in 
the right-hand sides of the constitutive laws obtained by dimensional reduction, 
see equations (63b). In particular, when the thread undergoes a rigid-body 



motion, uj_ is independent of S, the strain rate vector e vanishes and so does the 
viscous stress, as expected. 

Let us introduce the projections of e(5, t) in the tangential and normal di- 
rections, denoted et and gi, respectively: 



e{S, t) = e,[S, t) t_iS, t) + 6^(5, i), e^{S, t) ■ t{S, t) = 0. 



(35) 



The subscripts are motivated by the fact that the projections are associated 
with the twisting and bending deformations, respectively, as we show below. 
The projections are defined by the following explicit formulas, 



et{S,t)=e{S,t)-t_{S,t) 
e^{S,t)=P^{t_{S,t),e{S,t)) 



(36a) 
(36b) 



The covariant derivative in the right-hand side of the compatibility con- 



dition (33) commutes with the projection by equation (11). As a result, the 



decomposition of the twist-curvature vector tt in equation ( 30 ) implies 



etiS,t) = 
ey,{S,t) = 



dT{S,t) 
dt 

dK{S, t) 
dt ■ 



(37a) 
(37b) 



Here et and e^ appear to be given by the time-derivative of the strain measures 
— the kinematical twist r and the curvature binormal K^. These derivatives are 
evaluated in a frame following the material (notice the presence of the covariant 
derivative in the second equation, and recall that the time derivative in the first 
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equation is actually a convective derivative as we use Lagrangian variables). 
Equations ( 37 1 confirms that et and measure the rates of change of the twist 
strain r and of the bending strain K^, respectively. For the sake of completeness, 
we recall the definition of the rate of strain associated with the stretching mode 
of deformation, given earlier in equation (|20[) : 



d{S,t) 



d£{S, t) 
dt ' 



(37c) 



The decomposition of e in equation ( 35 ) , and the interpretations of its pro 



jections in equations (37) raise some difficulties in the discrete setting. Indeed, 



the definition of e as the gradient of rotation, given in equation (34), cannot be 
extended to the discrete case, as we shall see. However it is possible to propose 
discrete equivalents to its tangent and normal projections et and e^^ separately. 

3. 7. A geometrical identity for the twist 's rate of strain 

We derive in this section a geometric identity which is central to the mechan- 
ics of thin elastic rods and viscous threads. It explains the coupling between the 
motion of the centerline x{S, t) and the the twist t(S, t). In a previous work |35) 
focusing on the case of elastic rods, this equation was used to obtain a natural 
discretization of the twist. A similar discretization strategy is followed here. 



We start by projecting the compatibility condition (33) along the tangent 
direction, using the product rule: 



et 



dT{S,t) 
dt 



= t 



as 



d{t ■ u) _ dt 
~1)S dS 



By equation (27), we can identify the quantity (t-w) appearing in the right-hand 



side as the angular twist velocity v. In the second term, the time derivative of 



the tangent is given by combining equations ( 28 ) and ( 24 ) as 



et 



dv{S,t) 
dS 



ilL X t) 



In the right-hand side, we can replace tt by its transverse projection inside 
the cross product tt x We then permute the mixed product and use w x < = i 
by equation ( 25 ) . Then we find 



et 



dv{S, t) 
dS 



KiS,t) 



dt{S, t) 
dt ' 



(38) 



We shall now comment a little on this equation which is important both at a 
practical and at a fundamental level. 



On a practical side, equation (38) makes it possible to use the centerline 



position x{S,t) and the spin velocity v{S,t) as the primary variables for pa- 
rameterizing the thread. When combined the definitions I£ = tx t' and t — x' , 
equation ( |38[ ) provides the value of the strain rate for the twisting mode, et = t, 
a quantity that is required in the constitutive law (63b). This centerline/spin 



parameterization has important benefits, as we shall argue later on. 
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More fundamentally, equation ( 38 ) explains the coupling between the center- 
line motion and the twisting motion of the material frame: the twisting moment 
is proportional to the rate of strain et, which not only depends on the rotational 
degree of freedom v (through the first term in the right-hand side) but also on 
the centerline motion (through the second term). The geometry underlying 



equation ( 38 1 has been discussed by several authors but has never been used as 
a starting point for setting up simulations of viscous threads. This equation can 
be seen as an incremental version of the Calugareanu- White-Fuller (CWF) the- 
orem [57l [58l HH [60l m] which defines the notion of writhing for a closed curve 
— for a short review on this theorem, see references [SH [S3] . The CWF theorem 
has become very well known in the context of supercoiled DNA [331 ISH ISS [SB] 
or polymers |B7], and has been used in other contexts such as the dynamics 
of elastic filaments in a viscous fluid [54] [68] or the mechanics of proteins [69] . 



The binormal curvature K^{S, t) in the right-hand side of equation ( 38 1 above is 



directly responsible for the geometrically non-linear term T x n of the equations 



of motion, see equation (71). There is no obvious way to discretize this term 



of the equations of motion, but a natural discretization will become apparent 



when this term is connected to the geometric context of equation ( 38 1 , a variant 



of which will appear in the discrete setting, see equation ( 100 1 below 



3.8. Centerline /spin representation 

The initial parameterization of the thread introduced in sections |3.2| and |3.4 
is based on the centerline position x{S,t) and the orthonormal frame {d^{S,t)) 
with 1 < i < 3. This parameterization is subjected to the kinematical constraint 
of compatibility expressed by equation ( |24[ ) . From a numerical viewpoint it is 
desirable to eliminate this constraint, and use a reduced set of variables instead. 
To this end we introduce the centerline/spin representation {x,v): the center- 
line is described using x{S, t) as earlier but the material frame is described in an 



incremental way using the spin velocity v{S,t) introduced in equation (26). By 
'incremental', we mean that we do not keep track of the absolute direction of the 
transverse material vectors and but only of the spin velocity v = di-d^- In- 
deed, in the case of isotropic cross-sections which we consider here, the absolute 
direction of the material frame can be eliminated from the equations of motion. 
Our centerline/ spin representation is inspired from the centerline/ anf^Ze repre- 
sentation introduced by Langer and Singer in the context of elastic rods |53] : 
while they define the orientation of the cross-section incrementally with respect 
to arc length S using a twist angle 9, we define it incrementally with respect 
to time t using a spinning velocity v. The benefit of our representation is that 
it makes the matrix governing the dynamics of the thread sparse, as explained 
later. The initial (x, representation uses 3 degrees of freedom, such as the 
Euler angles, for the orientation of the material frame, which are subjected to 
2 scalar constraints for the compatibility. By comparison, the centerline/spin 
representation uses a single degree of freedom v{S,t) for the description of the 
twisting degree of freedom, which is free of kinematical constraint. This has the 
important benefit of decreasing the number of degrees of freedom and removing 
constraints. 
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In this section and in the next ones, we expose the principle of the time- 
stepping algorithm in the center line/spin representation. This algorithm is used 
at a particular time t to derive positions and velocities at the next time step. 
In the absence of any ambiguity, we shall often omit the time argument t from 
now on. The equations of motion being second-order in time, we assume that 
the actual position, as well as the linear and angular velocities as prescribed, 

x{S)^x{S,t), u{S)^x{S,t), v{S)^v{S,t). (39) 

We show how the linear and angular accelerations 7(5") = ^(5*) and v{S) can be 
computed. This requires reconstructing a number of intermediate quantities. 
To begin with, we can readily compute the axial strain i{S), the unit tangent 



t and the binormal curvature K_ defined respectively by equations (15), (16) 



and (31 ) in terms of the centerline, without even using the velocities: 



e{S) = \x!{S)\ (40a) 

m - (40b) 

K{S) = t{S) X t'{S). (40c) 

3. 9. Reconstruction of strain rates from velocities 

We now proceed to the reconstruction of the strain rates d, et and e^ defined 



in equations (37). We introduce a reconstruction scheme valid for arbitrary 
velocities, denoted u and that are not necessarily equal to the velocities u 
and V of the actual motion. These velocities u and v will be called virtual. 
Working out the dependence of the strain rates on arbitrary (virtual) velocities 
will enable us to put the constitutive laws in a variational framework, using 
dissipation potentials. This approach will allow for a natural discretization of 
the constitutive laws. 

We start by reconstructing the time derivative of the tangent, i. By equa- 



tion (21 1 this involves the operator V 



'S) = ^ Pi_ {m,u'iS)) . (41a) 

Note that V depends on the functions x and u, and not just on their values at 
S. In the right-hand side, £{S) and t{S) are reconstructed from the centerline 



x{S) passed as the first argument of V using equations ( 40 ) 



In the particular case of the real motion, when the operator is evaluated with 



the actual velocity u{S) = u{S) — x{S,t) defined in equation (39 1, V yields by 
construction 

mu;S) = ^^. (41b) 

Before considering the strain rates, we also need to reconstruct the material 



rotation lj. For this purpose, equation (27) suggests introducing the operator 



W{x; u, v; S) = v{S) t{S) + t{S) x V{x; u; S). (42a) 
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Here again, the tangent t is imphcitly a function of the first argument x. This 
operator W yields the material rotation in the case of a real motion u = u: 



W{x\u,v-S)=ui{S,t). 



(42b) 



The strain rates d, and are reconstructed by means of three operators, 
noted £s, -C* and respectively. The definition of the operator associated with 
the stretching mode is motivated by equation ( 20 ) : 

Cs{x;u;S) =t{S) -uiS). (43) 

When evaluated with a real motion it yields the strain rate, 

Csix;u;S)^diS,t). (44) 

The twisting and bending strains are introduced through a decomposition 
equivalent to the one used in equation ( 35 ) : 



C'ix;u,v;S)t{S) + €}'{x;u,v;S), 



dyVix;u,v;S) 
dS 

where the right-hand sides are defined by the following projections: 

dW{x;u,v;S) 



C{x;u,v;S) ^t{S) 
C^ix;u,v;S) = (t{S) 



dS 

dyV{x;u,v;S) 
dS 



(45) 

(46a) 
(46b) 



Note that the total derivative of >V in equation ( 45 ) and ( 46 ) produce terms pro- 



portional to u (S) and v'{S) according to the definition of W in equation (42a 
As a result, the operators £* and also depend on the local values of u and 



For a real motion, W_ evaluates to uj_ by equation (42b I, and so the left- 
hand side of equation ( 45 1 evaluates to the strain rate vector e. Comparison 
of the decompositions (35) and (46) shows that, by construction, the operators 
operators £} and £^ yield in the case of a real motion 

C\x;u,v;S) ^ et{S,t) (47a) 
£}'{x;u,v;S)=e^{S,t). (47b) 



A more explicit form of the operator C defined in equation (46a) can be 
found by inserting the definition of W into equation ( 42a I . After some algebra, 
we find 

£' (x; w, v; S) ^ v'{S) + K{S) ■ V(a;; u; S), (48a) 

where is the curvature binormal computed from the centerline configuration x 
using equation (40c I. This is an extension of equation (38 1 to virtual velocities. 

All the operators V, W, Cs, introduced above depend finearly on the 

virtual velocities u and v. Their discrete equivalents shall allow us to calculate 
the viscous stress. 
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3.10. Dissipation potentials 

In the Lagrangian framework, internal viscous stress can be described by a 
Rayleigh potential, see for instance reference [70]. This potential, which plays 
a role similar to the potential energy for elastic rods, expresses the power dissi- 
pated by viscosity during a virtual motion prescribed by the velocities u{S) and 
v{S). This potential has three contributions, corresponding to the stretching, 
twisting and bending modes of deformation: 

T^{x;u,v) = 'Ds{x-u) +'Dt{x;u,v) +'Db{x;u,v). (49) 

The stretching contribution is proportional to the stretching modulus whose 
expression is due to Trouton [T3] : 

Dii)=3^Aie) = ^, (50) 

where /i is the fluid's dynamic viscosity. Here Do = 3fi Aq = 3/i(7ra§) is the 
value of the stretching modulus in reference configuration. The twist modulus 
C and the bending modulus B read [7TJ |S] 

C{e) = 2^^I{£) = ^, B{£) = 3^,m = ^, (51) 

where Cq — 2/i Iq and Bq ~ 3 /i /q are the moduli in reference configurations, 
and the moment of inertia in reference configuration Iq is defined in section [33} 
We propose the following expressions for the stretching, twisting and bending 
contributions to the Rayleigh potential: 



Vs{x;u)^J^ ^^(^C,{x;u;S)ydS, (52a) 
V,{x;u,v)^ ^^^^(C\x_-U,v-S)f dS, (52b) 



V^{x-u,v)^J^ ?i^(^C\x;u,v;S)ydS. (52c) 

Here S" and 5*+ denotes the Lagrangian coordinates of the endpoints of the 
thread. Both 5*^ and 5^ may depend on time even though this time dependence 
is implicit for the sake of readability. In all the expressions above, £ is a function 



of the first argument x by equation ( 40 1 . Note that the stretching contribution 



Vs does not depend on the rotational degree of freedom v but solely on the cen- 
terline velocity u. Since Cs, C} and £^ are linear forms, all contributions I?s, 
and I?t and the total Rayleigh potential T) are quadratic forms of their velocity 
arguments u and v. This quadratic dependence reflects the linear character of 
the viscous constitutive laws. 



The expressions introduced in equations (52) for I?s, I?t and I?b will be 



justified a posteriori in section |4.2[ by checking that they lead to the correct 
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constitutive laws. In particular, we will explain the reason for the factor 1/t in 
the integrands above. 

The Rayleigh potential T) allows the equations of motion for a thin viscous 
thread to be put in variational form. As noticed by Batty and Bridson |72j in 
the context of 3D fluids with free boundaries, this variational setting provides a 
natural discretization of these equations. We follow this general approach, and 
start by exposing the variational structure of the smooth equations. 

3.11. Equations of motion 

The main property of the Rayleigh potential is that is gives the viscous force 
by derivation with respect to the virtual velocity: the resultant of the internal 
viscous stress on the centerline is given by 



OX>(x; M, v) 
du{S) 



(53a) 



The notation in the right-hand side must be understood as follows: we first take 
the functional derivative of the potential with respect to its argument u, and 
later substitute the velocity arguments with their real values, u = u and v = v. 
This quantity is the resultant of the internal viscous forces, per unit length 
dS in reference configuration. It includes the stretching, twisting and bending 



forces, each contribution being listed in equation (49). 

The stress conjugated to the spin velocity v is the twisting moment due to 
the internal viscous stress in the thread. It is given by a similar formula. 



QAs,t) = - 



dv{S) 



(53b) 

(u.v) — {u,v) 



In equations ( 53 ) the functional derivatives in the right-hand sides are calcu- 
lated practically by computing the first variation dV of the Rayleigh potential 
for small increments of the virtual velocities, denoted 5u and Sv, and by identi- 
fying the result with 

-dV{x;u,v;Su,dv) = (P^ Su + Sv) dS. (54) 

Js- 

This calculation will be carried out later in section [4^ when we work out explicit 
expressions for the net viscous force P.^ and moment Q^, and check that they 
are equivalent with the classical Kirchhoff equations. 

In terms of the net viscous force P^ and twisting moment Q^, the balance 
of linear and angular momentum can be written 

pAox{S,t)^P^{S,t)+P{S,t) (55a) 
£Jv{S,t) = Q^{S,t) + Q{S,t). (55b) 

Here P{S,t) is the density of external force and Q{S,t) the density of exter- 
nal twisting moment. These balance laws are for an infinitesimal segment, and 
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per unit length dS in reference configuration. The factors (pAo) and {£J) are 
its mass and moment of inertia about the tangent, measured per unit refer- 
ence length dS. Here, J is the moment of inertia per unit length ds in actual 
configuration, given by the usual formula 

J{£)^ II r^rdrde^2pl{£). (56) 

J J\7-\<a 



The factor £ in the left-hand side of equation (55b) is because an element of 



reference length dS* has a moment of inertia J ds — £ J dS. 



By multiplying both terms of the first equation (55a) by a vector- valued test 
function Su{S) and both terms of the second equation ( 55b[ ) by a scalar- valued 
test function 6v{S), we can write the equations of motion (55) in weak form. 
The motion is such that, for any choice of the functions 6u(S) and 6v{S) and 
at any time t, 

S+ y.S+ 

{pAou-du+Jvdv)dS = -dT>{x-u,v;du,Sv)+ {P-Su+QSv)dS (57) 



where u{S, t) = i{S, t) is the actual velocity. In continuum mechanics, the 
left-hand side is called the virtual work of acceleration, the first term in the 
right-hand side is the internal virtual work, and the last term is the external 
virtual work. This weak form of the equations of motion will be useful for going 
to the discrete case. 

3.12. External loading 



In equations (55), P and Q denote the the force resultant and the twisting 
moment due to external forces — as opposed to the internal, viscous forces in the 
thread. Those forces are given per unit length dS in reference configuration. In 
our validation examples, we consider a thread moving under the action of gravity 
and surface tension. Contact with the ground is not handled by applying forces 



but instead by freezing the motion of particles, as explained later in section 7.2 
we do not need an explicit expression for these contact forces. 
Gravity is represented by a force 

P{S,t)^P^{S,t)^pA„g, Q{S,t)^Qg{S,t)^0. (58) 

We now work out the expression of the effective forces acting on the centerline 
as a result of surface tension on the lateral boundaries of the thread. Those forces 
are derived by variation from an energy proportional to the lateral area of these 
boundaries. In the case of a slender thread which we consider, the capillary 
energy is written using the following approximation of the lateral area: 

£^{x)^ / 'y2TTa{£{S))£{S) dS, (59) 
Js- 

where 7 is the surface tension, possibly depending on time and position along 
centerline, and 2TTa{£) {£dS) is the lateral area of a cylinder of radius a and 
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length ds — £dS. We neglect the small conical angle of the lateral surface, 
which has a negligible influence on the capillary forces for a thin thread. 

The effective capillary force acting on the centerline can be obtained by 



variation, using the definition of a{£) in equation (23 1 and the definition of £ in 



terms of x{S) in equation (151. This yields 
d£{x; Sx) 



Ik-y Sx 



s+ 

s- 



PJS)-SxdS, 



(60a) 



where the bracket denotes the boundary term coming from the integration by 
parts, and 



dnJS,t) 



dS 



n (5, i) = 7 TT a{S, t) t{S, t). 



(60b) 

(60c) 
(60d) 



3.13. Neglecting rotational inertia 

For thin elastic rods or viscous threads, a classical approximation, proposed 
by Kirchhoff himself, is to neglect the rotational inertia, that is to set J = 
in the equation of motion (55b I. This approximation can be justified by 



the fact that the kinetic energy associated with rotational inertia scales like 
{£J)v'^ ^ £ pa"^ {l/t*Y for a motion happening on a typical time-scale t* . By 
contrast the kinetic energy associated with translation of the centerline scales 
like £ pAu^ ^ £ pa^ {L/t*)'^ , where L is the typical length-scale of the motion. 
The energy of the rotational mode is therefore negligible for slender threads, for 
which L 3> a. Rotational inertia is always dominated by translational inertia 
— except for a straight viscous thread moving in pure twist, a particular case 
where the kinetic energy in translation is exactly zero and the above argument 
is inapplicable. Therefore, we set J = in the following. The longitudinal 



balance of angular momentum ( 55b ) becomes a condition for the quasi-static 
equilibrium of the twisting mode: 



= Qv(5,t) + Q(5,i). 



(61) 



This equation expresses the fact that the typical time associated with the damp- 
ing of twist waves, which is much shorter than that associated with the damping 
of bending waves by the above scaling argument, is considered to be shorter than 
the time step of the simulation. 



4. Equivalence with Kirchhoff equations for a thin viscous thread 

This section aims at demonstrating the equivalence of the Lagrangian de- 
scription of threads based on the centerline/spin representation (x, w) exposed 
in the previous section, and the classical Kirchhoff equations for thin viscous 
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threads. The goal is to bridge the gap with classical formulations, and to iden- 
tify important stress variables that underlie the dynamic equations derived in 
the previous section. No new ingredient required in the numerical model will 
be introduced, and the reader interested only in the implementation can skip 
ahead to the derivation of the discrete model in section [5l 



4-.1. Constitutive laws underlying the dissipation potential 

To establish the connection between our formalism and the Kirchhoff equa- 
tions, we shall start by calculating the force and twisting moment Qv arising 
from viscous stress. By equation (54), this requires working out the first varia- 
tion of the dissipation potential with respect to the velocities u and v, near the 
real motion u — u and v = v. Combining equations (49) and (52), we can write 
this first variation as 



dl?(2;; u, v; 6u, Sv) 



£ 



dC^iSu^Sv) dS, 



where we have temporarily omitted the arguments x, u, v and S in the integrand 
for better readability. Since the latter are the real velocities, we can make use 
of equations (44) and (47), and write 



d'D{x; u, v; Su, Sv) 



Dd 



dCs{Su) 



—— dC {6u, Sv) + —— 



dC^iSu^Sv)] dS, (62) 



where the omitted arguments x, u and v again refer to the real motion. 

Let us denote rig the first coefficient appearing in the integrand, and assemble 
the two other coefficients into a vector denoted to: 



ns{S,t) 
rn{S, t) 



I 



{Dd) 

(Cett + Be^,) 



(63a) 
(63b) 



Note that these coefficients do not depend on the virtual motion but only on the 
actual one. The quantities introduced in equation (63) will be identified as the 
viscous stress in the thread. More accurately, is the scalar tension resisting 
stretching and to the internal moment arising from twisting and bending. 
The vector tension will also be useful later on, 



n^{S,t)^n,{S,t)t_{S,t). 



(64) 



It can be interpreted as the internal force that arises in response to stretching 
deformations. 
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4-2. Equivalence with Eulerian constitutive laws 

The constitutive laws can be rewritten in Eulerian form, which is how they 



are classically presented in the literature. From equations (63), we have 



n^{S,t) = D(f- 

rn{S, t) ^ [C {t(» t) + B {I - t(» t)] 



where we have introduced the Eulerian strain rate. 



1 

~l dS 



dui 



(65a) 
(65b) 



(65c) 



In the square brackets of equation ( 65b ) , we have introduced the tensor of 
viscous moduli. The operators in parentheses inside these square brackets are 
the tangential and perpendicular projection operators. 

In |Appcndix A[ we show that the constitutive laws used in the classical work 
of Ribe [9 , and derived by him from the Stokes equations in 3D, are equivalent 



to equations ( 65 1 above 



We note that, according to the Rayleigh- Taylor analogy, the case of an elastic 
rod is described by very similar equations, namely by replacing the strain rate 
in the constitutive law (65b) by the Eulerian twist-curvature vector tt^, and the 



constitutive law ( 65a ) for the internal tension by the condition of inextensibility, 
i=l. 



4-3. Canonical form of the internal virtual work 

The viscous introduced stress introduced in equations (63) allows equa- 
tion ( 62 ) to be written as 



— dV{x] u, v; 5u, Sv) = — / (^Ug dCs{6u) + m ■ dM{Su, dv)^ 



dS, 



(66a) 



where 

dM{x; Su, 5v; S) = t{S) dC\x; Su, 6v; S) + dC^ix; 5u, 5v; S). (66b) 

In this definition of the operator dAl, we recognize the first variation of the 
right-hand side of equation ( 45 ) — this can be shown by noting that its right- 
hand side is linear with respect to the virtual motion 6u and Sv. Inserting into 
the above equation, we find a compact expression for the first variation of the 
dissipation potential: 



dV{du, 6v) 



d[6u] 
"dS" 



m{S) 



d[dyV{6u,5v)] 
dS 



dS. (67) 



Here, the first term in the integrand has been rewritten using Ug dCg = (t 
du'{S)) = ■ Su'{S) by equations ^ and (|64). 
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For any real motion, the quantity W evaluates to the angular velocity lo 
of a cross-section by equation (42a I. Therefore dW can be interpreted as the 
virtual rotation associated with the virtual motion 5u and 6v. In view of this, 



equation ( 67 ) is similar to the classical expression for the internal virtual work 
in a 3D continuum, — /// g_: V(fe) d^x. In the context of thin viscous threads, 
the 3D stress tr is replaced by the stretching force in the first term and by the 
bending and twisting moment m in the second term, the gradient is replaced 
by the spatial derivative d/dS", and the virtual displacement {5x) is replaced by 
the virtual velocity 5u and virtual spin velocity 5v. 

4.4- Identification of the net viscous force and twisting moment 

We shall now derive explicit expressions for the effective force and moment 
acting on the centerline as a result of the internal viscous stress. Inserting the 



definition of V in equation (|41a| into the definition of W in equation (42a) and 

t X Su'{S) 



computing the first variation, we have 

dW{5u,Sv) =tSv + 



Inserting into the expression ( 67 1 for the virtual internal work, the first variation 

d{t6v+ ^tx Su'(S)y 



of the dissipation potential reads 



dV{6u, Sv) = 



dSu 



dS. 



ig_ \ dS dS 
We can integrate by parts to cast the right-hand side into a form similar to 



equation (54) 



- dV{Su, Sv) = BT 



S- 



dm , din^ + \txm'{S)) 

^'95^" + as 



where BT stands for boundary terms. Those boundary terms are omitted in 
the smooth case, and will be readily obtained from our variational formulation 
in the discrete case. 



Identifying equations ( 54 ) and ( 68 ) , we find an explicit expression for the 



net force acting on the centerline and for the net twisting moment 

dn{S, t) 



dS 



,iS,t) = t_iS,t) 



dm{S, t) 
dS 



where 



n{S, t) =n^ + 



l{S,t) 



t{S, t) X 



dm{S, t) 
dS ' 



(69a) 
(69b) 

(69c) 
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This net force results from the combination of viscous stretching, twisting and 
bending stresses inside the thread, and has been computed using the Rayleigh 
dissipation potentials. In these expressions we have omitted pointwise forces 
and moments at the endpoints S = coming from the boundary terms. They 
will be restored in the discrete model. 



4-5. Equivalence with Kirchhoff equations 

We proceed to show that the expressions (69) for the net force and 
twisting moment Qv acting on the centerline are equivalent to the equation of 
motion for a thin thread due to Kirchhoff. 

The Kirchhoff equations are usually written in Eulerian variables as 

dn 
ds 



P^ = pAx 



dm 



+ txn + Q^t = Jvt 



(70a) 
(70b) 



Here, P^ and are density of applied force and twist per unit length ds in 
actual configuration, respectively. 

When these loads are multiplied by ^ = we obtain our Lagrangian den- 



sities of load P 
by £, we have 



IP!^, Q = IQ^. Multiplying both equations (70a) and (70b I 



dn 
dS 



P_ = pAo'x 



dS 



+ Tx n + Qt = IJvt 



(71a) 
(71b) 



where Aq — £ A is the area of the cross-section in reference configuration, and 
T = It is the deformed (non-unit) material tangent defined in equation (14 1. 

Projecting equation (71b I along the tangent and normal directions succes- 
sively, we have 

dm 
'dS 



t ■ 



= £Ji 



P^{t,n) = \tx 



dm 
'dS 



(71c) 
(71d) 



The tangential component of the internal force n can be interpreted as that 
resisting stretching of the centerline. It is called the tension and is denoted n^. 
The full internal force n can then be reconstructed by combining this tangential 
component 



rigt with its normal component, given by equation (71d) 
1 dm 



(72a) 



Comparison of equations (71a) and (55a) reveals that the viscous stress, de- 
scribed by the quantities n and to, produces an net force 

dn 
dS 



P., = 



(72b) 
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Figure 3: Discrete setting: centerline is a polygonal curve. Note that we use subscripts for 
vertex-based quantities, such as vertex positions, and superscripts for segment-based quanti- 
ties, such as segment length 



on the centerline. Similarly, comparison of equations (55b) and (71c) reveals 
that the net twisting moment due to viscous forces reads 



dm 
~dS 



(72c) 



The expressions (72b) and (72c) for and Qv derived from the Kirchhoff 



equations ( 70 ) are identical to those derived earlier in equations ( 69 ) from our 
dissipation potentials, based on the centerline/twist representation. We have 
therefore established the equivalence of our formalism with the classical, Eule- 
rian equations for thin viscous threads. The benefit of our formalism is that it 
can be discretized in a natural way, and leads to an efficient implementation. 



5. Space discretization: the discrete viscous thread model 

In this section, the spatial discretization is carried out, by closely following 
the smooth formalism of the previous sections. Our derivation of the discrete 
model makes use of three key ideas that have been exposed in the smooth 
setting. First, we extend the centerline/spin representation to the discrete case; 
its benefit is to eliminate two out of three rotational degrees of freedom using the 
condition of compatibility of the tangent. Second, we introduce a discrete twist 
using the geometrical notion of parallel transport. Third, we derive equations 
of motion in the discrete setting by variational principles, starting from discrete 
dissipation potentials. 

We start our analysis of spatial discretization by defining discrete quantities 
such as centerline position, linear and angular velocities, rate of strain, etc. 
Time discretization will not be discussed until section |6l 



5.1. Kinematics of centerline 

The centerline is discretized using (n -I- 2) vertices. Their positions in space 
are noted XQ{t), Xi{t), . . . , a;„_|_i(i), see figure |3] Our numerical model involves 
setting up a force, assigning a mass, and integrating the fundamental law of 
dynamics at each vertex (t) . The thin thread behavior is produced by means 
of a discrete viscous force, which by design converges to the force P_^{S, t) defined 
in equations (53a) and (69) in the smooth limit, n — > oo. 
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The segment joining vertices Xj and x^^^i is noted 

r{t)^x,+i{t)-xM, (73) 

as shown in figure [3] FoUowing classical conventions, we use subscripts for 
indices < i < n + 1 associated with vertices, and superscripts for indices 
< i < n associated with segments. Since the vertex index i plays the role 
of the Lagrangian coordinate S, the segment vector T^{t) defined above is the 
discrete equivalent of the material, non-unit tangent T{S, t) defined in equa- 



tion (14). More accurately it is, like many other discrete quantities introduced 
next, an integrated quantity: the discrete tangent is essentially the smooth tan- 
gent multiplied by the discretization length. 

The discrete segment length P{t) and unit tangent i^{t) are defined by for- 



mulas similar to equations (15 T6| 

t{t) = \r{t)i (74) 

m = ^. (75) 

We define the vertex velocities by 

u_.{t) = ^M). (76) 

In terms of the velocities Uj{t), we define the integrated axial strain rate for 
segment T* : 

di'* (t) 

d\t) = ^ - t{t) ■ {u,+i{t) - uM), (77) 



in analogy with equations (17) and ( |20[ ). This is again an integrated form of 
the smooth strain rate d{S, t). 

The time derivative of the unit tangent is given in terms of the vertex ve- 



locities by a geometrical formula analogous to equation (21 ), 



W) ^Y(t)-^ {f{t),u,+,{t)~uM) ■ (78) 

To define the bending strain, we shall later need vertex-based tangents. 
There are several possible definitions that are equivalent in the smooth limit, 
and we opt for one that preserves the unit character of the tangent, namely 

Uix,_i,x„x,^:^) ^ ■ (79) 

The tilde notation is used here and in several other places when we introduce 
vertex-based versions of quantities that are primarily defined at segments, and 
vice-versa. 
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Similarly, there are several possible definitions for the discrete binormal cur- 
vature vector. One particular definition is; 



KM = , - ^ "^f (80) 
The motivation for choosing this particular definition comes from the forth- 



coming equation (100): this particular expression K^^ will emerge from the 



calculation of the discrete twist. The vector K^^ is an integrated measure of 



the smooth binormal curvature vector K_{S^t) defined in equation (311. In- 
deed the denominator in equation ( |80[ ) converges to 1 in the smooth limit 
where t*~^ ^ t' ~ t{S,t), while the numerator is equivalent to i*^^ x f 
X (f — f'~^) ^ K{S, t) £i where £i is the length of the Voronoi cell around 



vertex x^, defined below in equation (82 1 



5.2. Incompressibility: radius and related quantities 

Each segment carries a volume of fiuid and a mass of fluid m* . Those 
quantities are initialized based on the prescribed initial segment length, radius 
and mass density of the fiuid. They are conserved during the simulation, except 
in the case of an adaptive mesh, discussed in section [63] which requires segment 
subdivision. As in the smooth case we use incompressibility to reconstruct the 
local radius a*(t) and cross-sectional area A'^{t), assuming that each segment 
has a cylindrical geometry: 

AM-^y • (81) 

We shall need later the length £i of the Voronoi region near a given vertex. 
For an interior vertex with 1 < i < n, it is defined as the curvilinear distance 
between the midpoints of the adjacent segments, measured along the polygonal 
line traced out by the vertices: 

h{t) = ^''H^H^H^) forl<z<n. (82) 

This is a vertex-based discretization length, as opposed to the original segment- 
based discretization length £^ . This length is not required for the end vertices, 
i — Q and n -I- 1 . 

5.3. Material frame, angular velocity 

The unit tangent V; is defined at segments. We define the orthonormal triads 
(c?5^ , ^2 J ) a-t the segments too. This allows the condition of compatibility in 



equation (24) to be easily extended to the discrete case: 

m^t{t). (83) 
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Repeating the argument of section 3.4 one can show that the angular rotation 



oj* of the material frame can be decomposed as 

(i) = (t) f (t) + f (t) X f (t) . (84) 

In the first term of the right-hand side, the quantity v^{t) is the spin velocity, 
i. e. the angular velocity of the material frame about the tangent, as shown in 
figurejs] The second term warrants that the time evolution of the centerline, i 



0/ X remains consistent with the condition of compatibility in equation ( 83 ) 



5.4- Revisiting the case of zero twist: parallel transport 

The main difficulty in setting up a discrete model for thin viscous thread 
resides in the definition of twist. In the smooth case, twist is defined by project- 
ing the infinitesimal rotation vector tt along the tangent. This operation is no 
longer possible in the discrete case, as rotations are finite and are represented 
by a matrix. To remedy this difficulty, we introduce the geometrical notion of 
parallel transport. It enables us to revisit the notion of twist, in a way that 
makes its extension to the discrete setting natural. 

For a given configuration of the centerline, parallel transport defines a series 
of rotations from one segment to the next. Those rotations define a minimalist 
motion along the centerline, and will be used to define twist-less states of the 
thread. Parallel transport has also been used in the smooth setting to define 
the so-called natural or Bishop frame [SH |S3j . 

Consider the unit tangents and t^ of the segments adjacent to a vertex 
x^. We shall assume that these tangents are not opposite to each other, 

f-i ^ -t\ (85) 

This assumption is satisfied, except for a subset of configurations whose measure 
is zero. 

For a reason that will be clear in the next section, we define a rotation Q to 
be compatible at vertex if it maps t^^^ to t^: 

Q-f"i=f. (86) 

Next, we define parallel transport across vertex as the minimal rotation that 
is compatible; here, the word 'minimal' refers to the rotation having the smallest 
possible angle of rotation about its own axis Q This defines a unique rotation 



under the assumption of equation (85), as we show now. 

Parallel transport, defined above in geometrical terms, has an explicit repre- 
sentation: it is the rotation T. whose axis is along the binormal and whose 
angle is the turning angle (pi across vertex x.^. The turning angle is defined by 

(fii = cos^\f-'^ ■ f) with 0<(pi<TT. (87) 



^In geometrical terms, this minimal rotation minimizes the distance to the identity over 
the Lie group of direct rotations in Euclidean 3D space. 
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Note that (^i ^ tt by equation (85). The rotation just defined satisfies 

TJ ■T^=l (88a) 
T^-K.,=K, (88b) 
Z^-f'^^f (88c) 

These equations express the fact that T. is a rotation, that its axis is ahgned 



with binormal K^^, and that is it compatible — compare with equation (86). 
Compatibihty fohows from prescribing the rotation angle to be the turning 
angle ifi. In the particular case ipi = 0, that is when the adjacent segments are 
aligned, t^^^ = i', the three equations above no longer define a unique rotation; 
in this case, parallel transport is defined to be the identity, 

= 1 if (ySj = 0. (88d) 

Since a compatible rotation maps t^^^ and f, its angle of rotation has to 
be greater or equal to the angle (pi between them. The angle of rotation of the 
matrix T. that we have just defined is precisely ipi. Therefore, to show that 
the geometric definition parallel transport uniquely defines the matrix T , it 
remains to prove that any other compatible rotation has an angle of rotation 
strictly larger than ipi. This is what we do now. 

First note that any compatible rotation Q can be decomposed as 

g = Z^-E{f-\T{g)) (89) 

for some angle t{Q). Here i?(t*"^, r) denotes the rotation about t^^^ with angle 

T. This decomposition follows from the remark that T^^ • Q is a rotation leaving 

t^^^ invariant. Denoting q the unit vector obtained by rotating about the 
binormal by an angle cr, one can compute the dot product of q with its image 



• q by the rotation Q in equation ( 89 ) as 



q^ ■ q'^ — cos(pi + [am{ipi — cr) sin(o-)] (1 — cost{Q)). (90) 

This equality can be established in the direct orthonormal basis whose first and 
last vectors aret*~^ and K^/\K^\. respectively; in this frame, q^ = {cos cr, sin cr, 0} 
and — {cos(pi,smipi,0}. Details of the calculation are left to the reader. 
For any value of ifi such that < (pi < tt, there exists at least a value of cr 



that makes the function inside the square bracket of equation ( 90 ) negative and 
non-zero. If cost{Q) ^ 1, this implies 

q ■ q' < cos (fii. (91) 



If the rotation Q of equation (89) is different from T., cos t{Q) ^ 1. Then 



equation (91 ) shows that the angle of rotation of Q about its own axis is greater 
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or equal than cos ^{q^ ' '^l-)' ^'^^ therefore strictly larger than (pi. The proof 
is complete: there is a unique compatible rotation whose angle of rotation is 



minimal, and this is the rotation defined in equation (88). It will be called 
parallel transport. 

Parallel transport establishes a natural mapping between cross-sections be- 
longing to neighboring segments, and is similar to the notion of a Levi-Civita 
connection in differential geometry, see e. g. the book by Wald [73]. This prop- 
erty will now be used to define the discrete twist. More accurately, we shall 
define twist-less configurations of the rod to be those obtained by parallel- 
transporting the material frames from one segment to the next, and will define 
discrete twist by difference with parallel transport. 

The identification of parallel-transport with twist-less configurations of the 
rod can be justified by examining compatible rotations in the smooth case. In- 
finitesimal rotation tt are compatible when they can be associated with material 



frames satisfying the compatibility condition in equation ( 24 ) . Such vectors tt 



are of the form tt = X -I- r f by equation ( 30 1 . For a given configuration of 
the centerline, the binormal curvature X is prescribed but r is a free function. 
Smooth parallel transport is defined, as in the discrete case, by minimizing the 
magnitude |7r| of the infinitesimal rotation, keeping the centerline fixed; this 
minimization yields t{S, t) = 0. This result confirms that parallel transport is 
associated with twist-less configurations of the rod in the smooth case. In view 
of this, it makes sense to use parallel transport to extend the notion of twist to 
the discrete setting. 

5.5. Discrete twist 

Since material frames are orthonormal, there exists a unique rotation map- 
ping one material frame to the next. The rotation connecting the two material 
frames adjacent to the vertex x.^ is denoted Q , 

d)^Q^-d]r^ for j = 1,2, 3. (92) 
This finite rotation is the discrete equivalent of the infinitesimal rotation vector 



n{S,t) defined in equation (28); in the smooth case, the kinematical twist has 
been extracted from n by projection along the local tangent direction. This 
operation is no longer possible with the finite rotation matrix Q . 

Parallel transport provides an alternative route for defining twist in the 
discrete case. To begin with, note that setting j = 3 in equation ([92|) shows 



that Q is compatible — compare with equation (86), using cf^ — P. Being a 
compatible rotation, Q can be decomposed using parallel transport, as earlier 

g^ = 21^■ ^0 - i(f . n) ■ Z^- (93) 

Here, we use the notation — t{Q ). The alternative decomposition after the 



in equation ( 89 ) 



second equal sign in equation (93) follows from the following argument: the 
rotation [T. • R{f^^,Ti) ■ T^] is a rotation of angle t^, being conjugated with 
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R{f ^,Ti), and in addition leaves invariant — it is therefore the rotation 
about f' with angle r^, which is noted t^). Note that the angle Ti is uniquely 
defined modulo 2tt. 

In equation (93), the angle defines an axial rotation required to match one 
material frame to the next, in complement with parallel transport This Ti 
is called the discrete angle of twist across vertex Xj. It is an integrated version 
of the kinematical twist t(S, t) appearing in the smooth setting. 

Let us check that the discrete twist Ti is consistent with the smooth twist 
T{S,t) in the smooth limit. In equation (93), the rotations Q , T. and R{V-,Ti) 

|_ =i — * — 

converge in the smooth limit towards infinitesimal rotations which are repre- 
sented by the vectors n, K_ and Ti t respectively; the proof of this is left to the 
reader. In this limit, equation (93 1 becomes tt = + r^i, and we recover the 
decomposition in equation ( 30 ) . This confirms that our definition of the discrete 
twist is consistent in the smooth limit. 



5. 6. Rate of change of twisting strain 

To simulate the dynamics of viscous threads, we need an expression for the 
rate of strain associated with the twist mode. By analogy with the smooth case, 
it is defined as the material derivative of the angle of twist, 

e\ = h. (94) 

Like the quantity r.;, this is a spatially integrated form of the smooth rate of 
strain et(S', t). The goal of the present section is compute e\ in a form that can 
be used in our centerline/spin representation, i. e. to express e* as a function of 
the velocities and . 

To this end, let us start by introducing the polar angles r~ and t^ of the 
binormal K_i in the frames (di~^,d2^^) and (^5^,(^2)1 respectively. These angles 
are represented in figure |4] They can be computed by taking the arctangent of 
the coordinates of the vector t!^~^ x t*, which is aligned with K_i by definition 
of the latter. These coordinates are denoted {t^^t^^) in the first frame and 
(Tj^, t^) in the second frame: 

T± = tan-1 (^^^ , (95a) 

where 

r± = (f-i X f ) . (95b) 

In the second equation, j takes on the values 1 or 2, and the superscript in the 
last factor evaluates to the index (« — 1) of the segment in the left-hand side 
when ± = — , and to the index (i) of the segment in the right-hand side when 
± = +. 

As shown graphically in figure |4j the rotation Q can be decomposed into 

a rotation about f'^^ with angle t^ that brings (^['^ onto the binormal Kj, 
composed by the parallel transport that maps t^~^ to t^ without affecting the 
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Figure 4; Illustration of equation | |96[ l. The rotation mapping one material frame to the next 
is decomposed using parallel transport and the polar angles t~ and . 

binormal, and composed by a rotation about V; with angle {—t^) that brings 
back the binormal to d^i without affecting the tangent: 

Q=R{t\-rt)-T^-R{t-\T-). (96) 



As earlier in equation ( 93 1, we can use conjugation to group the axial rotations in 



equation above. Identifying the result with the definition of in equation (93) 
we have 



T. = T; - Ti 



Inserting into the definition ( 94 1 of the strain rate, we have 



(97) 



We proceed to compute the time derivatives by differentiating equations ( 95 ). 

The calculation of f~, is done in the frame moving with the first ma terial 
frame (c^}~^)i<j<3- There, all vectors in the right-hand side of equation (95b I 



are still, except which has angular velocity — . This yields, for j = 1, 2, 

= (f-i X [(w' - X f ]) • = -(f-i X d]r^) ■ [{u' - u'-^) X f ], 

after permutation of the mixed product. Inserting this expression into the 



derivative of the arctangent in equation (|95a|) , we find 

''il "^12 ~ '''■12 '''il 



r.a r^2 - ^.2 ^ (f x f ) ' [(^^ - ^''^) x f 



(98a) 



The time derivative of the second angle is given by the same formula, with 
the indices i and i ~ \ swapped: 

. + _ (f-ixf).[(c.'-c.'-i)xf-i] 
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Inserting this expression into equation (97) and permuting the mixed product, 
we have 



„t 



X f |2 



In the right-hand side, the second factor can be simphfied using the fact that 
both and f are unit vectors. This yields 

(99) 



Inserting now the decomposition (84) of the material velocity uj^ and sim- 
plifying, we have 



v'-'+Kr'^^^, (100) 



after using the definition (80) of the discrete binormal curvature K^^. 

This formula is fundamental as it yields the rate of strain for the twisting 
mode in a form that is suitable for our centerline/twist representation. It closely 



resembles the smooth equation (38). The second term in the right-hand side 
has a geometrical origin. It captures the change of parallel transport resulting 
from a change in the centerline, an effect that was dubbed holonomy in our 
previous work. The holonomy term is responsible for the coupling of centerline 
motion with the twisting mode, a phenomenon which appears to be geometrical 
in essence. 

The geometrical origin of this coupling has been recognized earlier but has 
not been used as a starting point for dynamical simulations of threads, to the 
best of our knowledge. The role of the binormal curvature K_ has been noted 
in the related context of Fiiller's theorem [61, for the increment of writhe of a 



space curve. Expressions similar to those in equations (98 1 have been derived 
for the increment of discrete writhe, and have been used for the simulation of 
the Brownian dynamics of DNA modelled as an elastic rod [74] . 

5. 7. Rate of change of bending strain 



In the smooth case, we have defined in equation (34) the strain rate vector 
e{S, t) to be the gradient of rotation. We introduce similarly a discrete strain 
rate vector e* by 

e'=ui'-ui'-\ (101) 



In equation (35), we have shown that the tangent and perpendicular projec- 
tions of e are the rates of strain relevant to the twisting and bending modes, 
respectively. In the discrete case, we carry out a similar decomposition, using 
the vertex-based tangent t^, 

el^h,{ip,)i-e il^K{v,)P^{l,e,). (102) 

In these projections, have introduced two normalizing functions ht and ft-b of 
the turning angle ipi. These functions reflect the fact that the definition of the 
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vertex-based tangent in equation ( 79 1 is somewhat arbitrary. For instance. 



we could use instead a different vertex-based tangent, 

1= ^^l.X^^ =H^i)k- (103) 

Here h{ipi) = ||^| = 1/ cos((pi/2), as can be shown by using trigonometric re- 
lations in the triangle whose sides are t^~^ and t^. For consistency with the 
smooth case, we require that the functions ht and /ib converge to one when 
their argument vanishes. 

Our definition of a discrete twist based on parallel transport imposes a par- 



ticular choice of the function h^. Indeed, identifying equations (99 1 and (1031 
we find e ■ = • e and so ht{(pi) — h{(pi) — ||^| — 1/ cos{ipi/2). 

By contrast, any choice of the function h\i{ipi) is acceptable as long as ft,b ^ 1 
for Lpi — > 0. This leads to infinitely many different discrete thread models, which 
are all equivalent in the smooth limit. In this paper, we choose /ib(</'i) = 1 for 
simplicity: the strain rate associated with the bending mode then reads 

i^ = PAke;). (104) 

5.8. Reconstruction of strain rates from velocities 

As in the smooth case, we keep track of the dependence of all secondary 
quantities on the velocities, and introduce virtual vertex velocities and vir- 
tual spin velocities at the segments. This shall enable us to compute the 
discrete viscous forces as the gradient of the dissipation potential with respect 
to velocities. 

The vertex positions are collected into a generalized coordinate Xit)^ a vector 
of dimension 3(n + 2): 

m^{^0{t),---,^n+l{t)}- (105) 

Note that there is no need to keep track of the orientation of the material frame 
in the generalized coordinate X(t), as we consider isotropic cross-sections. The 
twisting mode will only be included in the generalized velocity C/(i), defined 
later, through which it gets coupled with the centerline motion. 

Given the centerline configuration X(t), we first compute all quantities, such 
as P'{]C), ti'iK), iiiK) and Ki^{X_), which do not depend on velocities. To make 
the notations lighter, the argument 2L will often not appear explicitly in the 



following. Next, we extend the operators V and >V defined in equations (41 1 



and ( 42 ) to the discrete setting. The operator V' is attached to segment T' and 
defined by 

F(^;m„M,+i) = Y,Pi_ (f -M,) . (106a) 



This definition comes from equation (78): by design, this operator yields the 



time derivative of the tangent when applied to a real motion, as in the smooth 
case. 
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The discrete operator associated with the segment is defined by 



ii+l 



(106b) 



This definition is motivated by equation (84): when evaluated with a real mo- 
tion, W' yields the angular velocity vector u/. 

We now propose discrete versions of the three fundamental linear forms 
defining the viscous dissipation potential, which have been introduced previously 
in equations ( 43 1 and ( 46 1 . 

The discrete axial strain rate (T on segment T* is given by equation ( 77 1 . In 
view of this, we introduce the linear form 



(107a) 



Then, we have Cl{X_;u^,u^_^^) — t — for any real motion. 

The discrete strain rates for the twisting and bending modes are given by 
equations ( 100 1 and ( 104 1, respectively. Dependence of these strain rates on the 
velocities is captured by the following operators 



V — V 



(107b) 



and 



(107c) 



For any real motion, £* and €^ are equal to the strain rates e* — Ti and e^, 
respectively. The definition makes use of the geometrical definition of discrete 
twist based on parallel transport. 



5.9. Dissipation potentials 

In our centerline/spin representation, the generalized velocity is a vector of 
dimension 4n + 7 defined by collecting the linear velocities at the vertices, and 
the angular velocities of spin at the segments: 

U{t)^{uo{t),v°{t),u,{t),v\t),--- ,v"{t),u„^,{t)}. (108) 

This vector C/(i) is the generalized velocity of the real motion. The dissipation 
potentials are defined in terms of arbitrary velocities, which we collect under 
the name of generalized virtual velocity f/: 

mt) = {Mo(i),^^°(i),Mi(t),*'(i),--- ,i)"(i),M„+i(i)}. (109) 
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As in the smooth case, the viscous internal forces are introduced by means of 
dissipation potentials. The discrete potentials extend the smooth ones defined 



in equations (|52 1 

n 



0<i<n 
l<i<n 



(110a) 
(110b) 
(110c) 



i<j<ji 



Note that the stretching contribution involves a sum over all segments, although 
the twisting and bending contributions involve a sum over interior vertices. This 
is because the discrete axial strain is defined on segments, while the strain 
rate vector relevant to the twisting and bending modes is defined on interior 
vertices. 



In equations (110), the discrete moduli are defined by 



3fi* A' 

I. 



(111a) 
(111b) 

(111c) 



where is the fiuid's dynamic viscosity which is stored at segments like other 
fluid pro perties, is the segment's cross-sectional ar ea r econstructed by equa- 
tion (81), the segment length given by equation (74) and t the length of 



the Voronoi cell around an interior vertices given by equation ( 82 ) . The factor 



appearing the twisting and bending moduli is defined at vertices by linear 
interpolation over the adjacent segments: 



1 (A*-l)2 

2 



An 



(Hid) 



This definition is motivated by the fact that / = A'^ /{Att) in the smooth case, as 
shown by equation ( 22 ) . Note that the discrete moduli satisfy the same relation 



Bi/Ci — 3/2 as in the smooth case; physically, this relation is a consequence 
of the fiuid's incompressibility. All moduli depend on the actual configuration 
X{t) but not on velocities. 

The definitions (111) of the discrete moduli are identical to the defini- 



tions (50) and (51) in the smooth setting, up to factors proportional to the 
discretization length P or ii. These factors were introduced so as to warrant 
convergence of the dissipation potentials in the smooth limit. For instance, for 
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the stretching contribution we have ~ and CI ^ Cg t by equation ( 107a 
Formal convergence of the corresponding dissipation potential follows: 

i i i 

The total dissipation potential is defined by summing up the stretching, 
twisting and bending contributions: 



v{x-u) =v,{X:,u) +Vt{X:,u) +v^{x-u). 



(112) 



5.10. Discrete equations of motion 

The Rayleigh dissipation potential V is used to derive the discrete viscous 
forces and moments. By analogy with equations (53a I and (53b), we model 
the internal viscous stress by a net force acting on the vertex x^, and by a 
twisting moment Qt, acting on the segment Tj' , both of which are given by a 
derivative of the dissipation potential: 



dV{X]U) 



dV{X]U) 



u=u 



u=u 



(113a) 
(113b) 



In these expressions, the partial derivative is with respec t to the or entry 
inside the generalized coordinate vector [/, see equation (109 1. 
The discrete equations of motion read 



m.,u^{t) ^ p:{x{t)m)) + m) 

t r v\t) ^ Ql{X{t);U{t)) + Q\t), 



(114a) 
(114b) 



where P_i{t) and Q^{t) define the external loading, ihi is the vertex-based mass, 
defined as the sum of half the mass of the segments adjacent to vertex x^, 



E 



2 ' 



(115a) 



where the set of adjacent segments is Jj = {i — 1, z} if < z < n + 1, Jq = {0} 
and Jn+i = {'^}- Note that this vertex mass could change over time if 
adaptation were used. In equation (114b), is the density of moment of inertia 
of the cylinder attached to segment T* in actual configuration, per unit length: 



J' = 2«' r 

and — m^/V is the mass density of segment i. 



(115b) 
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The equations of motion can be written in compact form by introducing 
the generaUzed viscous force -F^, the generahzed external force F and the mass 
matrix M. The latter are obtained by collecting vertex and segment-based 
components with the same ordering convention as in the generalized velocity U: 



E-v — (£o tQvi E-lT ' ' ' Qv 1 £n+l ) 
Z = {Pq, Q",Pi, ■ ■ ■ Q",Pn+l) 

M = diag(TOo J°,mil,--- , T J 



" m„+il), 



(116a) 
(116b) 
(116c) 



where 1 represents the unit matrix in 3 dimensions. The equations of mo- 



tion (114) and (113) can be rewritten as follows, 



M-iLit) = F,{x{t).,u{t))+m 



m - n 



dU 

u{t). 



(117a) 
(117b) 
(117c) 



The third equation is the definition of vertex velocities, Ui{t) = ii(t). This 
equation makes use of the projection operator mapping the degrees of free- 
dom associated with vertices from the X representation, from which segments 
are absent, to the representation: 



ri+l 2 
i=0 j=0 



lii+i ■ 



(118) 



Here is a matrix of size (3ri -f 5) x (4n -I- 7) , defined with the convention that 
vector indices start at 0. The index i runs over vertices, the index j over space 
directions, and 5^ represents the vector whose entries are all 0, except for the 
fc-th entry whose value is 1. The values k — ^li -\- j and k = 4i + j appearing in 
subscript are the indices of the degrees of freedom for vertex i in direction j in 
either representation. 

Since the discrete dissipation potential is consistent in the smooth limit 
n — oo, the discrete viscous thread model converges to the smooth model in 



this limit: formally, the dynamical system in equations ( 117) becomes equivalent 



to the smooth equations of motion ( 55 ) combined with the expression n53« for 



the internal viscous forces. This convergence is checked numerically in section [Sj 



5.11. Surface tension and other forces 

The weight of the thread is taken in account by setting P,- — grhi and 
= in the equation of motion (116b). Here nii is the mass attached to a 
vertex, defined in equation (|115a), and g the acceleration of gravity. In our 
validation experiments, the thread is formed by expelling fluid from a syringe, 
and letting it fall onto the ground under its own weight. The contact forces 
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Figure 5: Surface tension is based on a cylindrical representation of the fluid attached to 
segments. 

with the syringe or the ground need not be computed, as we treat contact 
using kinematical constraint and not penalty forces. This is explained in the 
forthcoming section |7.1| 

We shall now explain how surface tension is taken into account — this is the 
last type of forces that we shall need in our examples. The present implemen- 
tation of surface tension assumes that the segments are cylinders, as in figure [Sj 
It is significantly simpler than that based on truncated cones presented in our 
conference paper |41| . The lateral area of the cylinder joining vertices and 
a^i+i reads 

= 27r a* f = 2 n/^Tv^, (119) 

as can be shown using — n {a"^)^ t. Its gradients with respect to the position 
of its endpoints read 

— ii+l 

Here we have used V^, = V^ lxj^^ — = — and V^. = +f. Discrete 
surface tension forces are set up by means of the discrete capillary energy, 

n 

£,m^Y.^''^\x). (121) 

Here 7' is the fluid's surface tension at segment i. In this equation, we assume 
that the radius a* of the cylinder varies over much longer length-scales that than 
the radius itself, and neglect the longitudinal curvature of the lateral boundary 
in front of its azimuthal curvature. This is consistent with the thin thread ap- 
proximation used everywhere in this paper. Let us note this approximation is 
not suited to the analysis of the Rayleigh- Taylor instability, whose critical wave- 
length is comparable to the radius. This instability would have to be studied 
using the full equations for 3D viscous fluids anyway, and not the dimensionally 
reduced equations for thin threads. 



1/2 



4-i „ ^1 J.2 

t = — TT a t 



1/2 



f = +TraH\ 



(120a) 
(120b) 
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The expression for the discrete capillary forces at a vertex is given by minus 



the gradient of the capillary energy (121) with respect to vertex positions 



-n^t f if i = 

P] = { if i = n + 1 (122a) 

nj^tf -n^i^f-i ifl <i<n 

Ql, = 0. (122b) 

The second equation is a consequence of the fact that the capillary energy 
depends only on vertex positions and not on the twist degree of freedom. 

These vertex forces P] are caused by an longitudinal force nl^^, called the 
line tension, acting along each segment. Its magnitude is given by identification 



with equation (120), 

n^,^7rjW. (123) 

This expression for the line tension nl^^ can be interpreted as resulting from the 
overpressure in the fluid caused by the interface curvature 1/a*, according to 
the Young-Laplace law. As it derives from energy proportional to the lateral 
area, the capillary force tends to make the thread shorter and more compact, by 
bringing the endpoints closer to each other and by flattening out curved regions 
of the centerline. 

The net capillary force on a vertex given in equation ( 122[ ) has two contribu- 



tions that almost cancel each other at each interior vertex (these contributions 
are associated with each one of the two adjacent segments), but only one contri- 
bution at the terminal vertices Xq and x^j^^. This reflects the presence of Dirac 



contributions at the endpoints in the smooth model, see equation (60a). Note 
that the area of the cap closing the cylindrical thread near its endpoints is neg- 
ligible, and the Dirac contributions are accurately captured even though these 
caps are not taken into account. Since the discrete model for surface tension has 



been derived from the energy in equation ( 121 ) which is a good approximation 



of the capillary energy in equation (59), the discrete capillary forces converge 



to the smooth ones in the limit n ^ oo. Our discrete model for surface tension 
is validated in section l8?3l 



6. Time discretization, numerical implementation 

6.1. Representation of the dissipation potential by a band matrix 

Th e str ain rate operators Cl{X', IL), ^li^'^ IL) and C!^{X.', LL) defined in equa- 



tions ( 107 1 are linear with respect to their virtual velocity argument IJ_. The 
two first operators are real-valued: each one can be represented as a vector, 
denoted ^s(^) or £-(X), that acts on [/ by dot product. The last operator is 
vector- valued, and can be represented as a matrix £}^{X) acting on U_: 

Cl{X;U)^Cm-U (124a) 
Cl{X;U)=£l{X)-U (124b) 

^■(^;m=^l'(^)-^- (124c) 
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The tensors Cl{X), ^liK) and €}^{X) just introduced can be distinguished 
from the original functions as they bear one additional bar below, and have a 
single argument. They can be calculated explicitly (?) by setting the virtual 
velocity t)_ = 6_j, i. e. by canceling all entries of U_ except for the entry at index 
j which is set to one, (ii) by evaluating the left-hand sides above using the 



definitions (107), (in) by filling the entry or column with index j in the target 
vector or matrix, and (iv) by iterating over the index j. In other words, we 
build the tensors using 



4n+6 
3=0 

and similar formulas for Cl{X) and C^{X)- 

Since the linear forms Cl{X_;l)_), C\{X;U_) and ^{X;U) attached to a seg- 
ment or a vertex i depend only on the virtual velocities of the neighboring 



segments or vertices, the vectors and the matrix introduced in equations (124 1 
are sparse. This sparse structure allows for efficient storage and manipulation. 

Similarly, the dissipation potential is a quadratic function of the virtual 
velocity V_. It can therefore be represented by a symmetric matrix VlX), such 
that for any virtual velocity U_ 

'DiX;!/) = ^U-2,{K) -LL- (125) 

This matrix is built up by combining the stretching, twisting and bending con- 
tributions, 

E(^) = a® + 2t® + 2b®- (126) 
Explicit formulae for these contributions can be found by inserting the repre- 



sentation of the linear operators in equations (124 1 into the definitions (110) of 
the discrete dissipation potentials: 

1.®= E D\X)C{X)^ cm (127a) 

0<i<n 

P(X)= C,{X)£}AX)®^l{X) (127b) 

l<i<n 

l<j<Tl 

Because the linear forms in the right-hand sides are represented by sparse ten- 
sors, these symmetric matrices in the left-hand sides are all band-diagonal. Their 
band structure is shown in figure |6] 

Note that new vertices are created during the simulation, while others collide 
with obstacles and are discarded. As a result, the number of vertices — and 
therefore the dimension of the dissipation matrix — may vary from one time 
step to the next. 
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Figure 6: Band structure of the dissipation matrices 2? {X_), 'D^{X_) and ^^(X} for the stretch- 
ing, twisting and bending modes. 



6.2. Boundary conditions and kinematical constraints 

We consider the possibility that some degrees of freedom are kinematicaUy 
constrained. Such constraints are typically used to enforce boundary conditions. 
For instance, in the case of a clamped end, the velocities Uq and of the first 
two vertices and the angular spinning velocity of the first segment v'^ are imposed 
by the motion of the clamp. 

In the presence of kinematical constraints, the the generalized velocity C/j^j 
at the end of the time step is prescribed to be of the form 

U,+,=R-W^+, + B', (128) 

where t is the time at the beginning of the time step, e is the time step, and W^^^ 
is a reduced velocity vector collecting independent degrees of freedom. This 
vector W_f^^ is the unknown of the time-stepping algorithm. Its size, denoted 
r, may vary from one time step to the next as constraints can be created or 
destroyed — this happens for instance when the thread makes its first contact 
with the ground. 

The matrix i? is a (4ri-t-7) x r matrix that dispatches the independent degrees 
of freedom into the generalized velocity vector U_. Given a strictly increasing 
numbering b of the unconstrained degrees of freedom, 6(0) < 6(1) < • • • < 
6(r — 1), this matrix reads 

r-l 

t=0 
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where is the Dirac vector, {5.j)i — 1 for i = j, and otherwise. 



The last term in equation ( 128 ) is the vector 13' collecting the velocities of the 



constrained degrees of freedom in the final state U_t+e- This vector is typically 
filled using the motion of the bodies in contact with the rod as explained in 
sections [Tj Entries of S' corresponding to unconstrained degrees of freedom can 
be set to zero by convention: 



■ B' ^ 0. 



(130) 



Indeed, allowing them to be non-zero would simply shift the value of the un- 
known 

As an illustration, consider the case when there is no kinematical constraints. 
Then r = An+7, b{i) = i and 13' = 0, which implies that B is the identity matrix. 
The more interesting case of a viscous thread, both ends of which are attached 
to clamps having a motion of pure translation with velocities ul\i^„^p and "clamp 
is handled by setting 



r = 4n + 7- 2x7 = 4n-7 
b{i) =i-7 

— ^ (—clamp' 0, Mclamp' 0' Qj ' ' ' i Q, 0, Mclamp' ^7 Mclamp)' 



(131a) 
(131b) 
(131c) 



In the presence of kinematical constraints, we discard the equations of motion 
corresponding to the constrained degrees of freedom. This is achieved by left- 
multiplying both sides of equation (117a) by B^ , and the equation of motion 
becomes 



^ ■M-U{t)^B^ ■ {F^{X{t),U{t)) + F{t) 



(132) 



6.3. Time discretization and time- stepping 

The time step e can be fixed or variable. At each time step, the updated 
position 2Lt+e aiid velocity f/^+g must be determined from the actual position 
X f and velocity U_^. We discretize equation (132) in time using a linear implicit 



scheme. The viscous force -F^ is evaluated implicitly with respect velocity but 
explicitly with respect to position: 



B'^ -M ■ 



U 



t+e 



= ■ {F,{X„U,^J + F{t)) ■ 



(133) 



This choice combines good stability, as demonstrated by the validation exam- 
ples, and ease of implementation: only a linear solver is required. 

The viscous force F^ (2Lt j ILt +t) i s indeed linear with respect to the unknown 
ILt+e by equations (117b) and (125), 



F^{x:,u) = ~ = -mx)-u. 



(134) 



Inserting equations (128) and (134) into the equation of motion (132), we find 



that the update rule for the velocity takes the form of a linear equation to be 
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solved for W. 



-t+e- 



■W 



t+e 



■ (6 i-mx,) ■g + F{t))-M-{B'-U,)). 

(135) 

In this equation, i? and B' encode kinematical constraints, M is the mass ma- 
trix, 2? the viscous dissipation matrix representing internal stress, F[{t) is the 
external loading and the velocity at the start of the time step. The matrix 
in square brackets in the left-hand side is symmetric, positive definite for any 
value of the time increment e > 0. As a result, this linear equation can be solved 
using efficient and robust solvers. 

Note that the external force -F(i), which includes in particular the effect of 
capillary forces, is evaluated explicitly. Wc have tried a linear implicit imple- 
mentation of surface tension, but have not observed any significant improvement 
in stability: an explicit evaluation of surface tension force does not raise any 
stability issue in all our tests. 

Once the linear equation for Wf^^ has been solved, the generalized velocity 
JZt+e is reconstructed by means of equation (128). Positions are then incre- 



mented using a discrete version of equation (117c 



n -u 



t+e- 



(136) 



6.4- Time- stepping: summary 



The main task of the time stepping loop is to set up and solve equation ( 135 ). 



Algorithm [T] explains how the various quantities appearing in this equation are 
constructed, and provide an overview of the implementation of the time step. 



Require: e {time step} 
Require: m*, V^, ji^ {fluid properties} 
Require: x^{t), Ui{t), v^{t) {initial positions and velociti es} 
Require: B, B' {kinematical constraints at boundaries, [6.2 {7^ 



Require: Frp {external force, j |5.11 



and eg. (116b)} 



md i 



5.21 



assemble 2Lt and jeqs. ( 105 ) a nd (|1U8 )} 



set rhi and J' {eqs. (1151 and/or (137)} 
assemble M {eq. (116c)} 



solve for W_f_^_^ {eq. (T35)} 



set i\ and £i {geometry, ^5.1 
set D^, d, Bi {viscous moduli, eq. ( 111 \} 
set £g, £^ {discrete strain rates, [5.8 
set V {dissipation matrix, eqs. ( 126 ■ Tiff} 



and S 



6.n 



reconstruct C/j+e {eq. (128)} 



update u^(t + e), u*(i +7){eq. ( |108| )} 
update Xj(t + e) {eq. (136)} 



Algorithm 1: Overview of a time step. 



Note that fluid properties such as the mass , volume and viscosity /i-' 
are stored in segments. These properties are constant in most of the examples 
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shown here. They may however depend on time, either when adaptation is used 



(see section 6.5), or when the dynamics of the thread is coupled to another 
physical process, as for instance in the case of heat transfer and temperature- 
dependent viscosity — a separate update rule is then required. In addition 
segments store their spinning velocity . 

By the argument of section |3.13[ rotational inertia can be neglected for a 
thin thread. In all the examples shown here, the rotational inertia was set to 
zero, 

= 0. (137) 
We have compared simulations done using this approximation with simulations 



based on the exact value of given in equation (115b), and checked that the 



results are indeed very close when the thread's radius is small. 

6.5. Adaptive mesh refinement 

In the the experiments of Morris [13], which we reproduce in section 9.2 
gravity stretches the thread by a factor which can be as large as 10 to 100. 
In the absence of refinement, the segments in the bottom part of the thread 
would be considerably longer than those at the top. A good spatial resolution 
is needed at the bottom, where the impact with the hard surface forms a coils 
of small radius. This makes simulations of severely stretched threads extremely 
inefficient, unless spatial refinement is used. 

We have implemented adaptive mesh refinement as follows. At the end of 
every dynamic step detailed in algorithm [TJ the segments needing refinement 
are first marked according to some user-defined criterion; next, segments that 
have been marked are actually split: a new vertex is inserted in each segment 
undergoing subdivision, and quantities such as mass, radius, position, velocity 
etc. are calculated in the new vertices and segments, as explained below. The 
refined thread is then used as the initial state for the next dynamic step. 

The refinement criterion can be based on comparing length of the segment 
to a maximum prescribed length. Other criteria based on the turning angle 
with respect to the neighboring segments have been considered too, but not 
used in the examples. The maximum allowed length may be a function of the 
position of the midpoint of the segment, to force refinement near a boundary for 
instance. In each case, the refinement criterion was adjusted manually to offer 
the best compromise between efficiency and accuracy — even if introduced in the 
context of non-steady problems such as the viscous sewing machine, refinement 
was always first validated in a steady coiling geometry, as explained in section [Sj 

Whenever a segment has been marked, its subdivision is carried out as fol- 
lows. A new vertex is inserted, two new segments are allocated and the former 
segment is removed. The position x.^ and velocity of the new vertex are cal- 
culated by an interpolation of order 4 based on the positions and velocities of 
its neighboring vertices. Vertices resulting from a subdivision concomitant with 
that under consideration do not enter in this interpolation, and so the result is 
independent of the order in which marked segments are processed. The mass 
stored in the original segment is equally split among the two subsegments. The 
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viscosity /i^ , surface tension 7^ , spinning velocity of the new subsegments are 
all computed by an interpolation of order 2 based on the values of the former 
segment and of its neighbors. Finally, the volume of the subsegments is 
computed by first considering an interpolation of the cross-sectional area at 
order 4, which is then multiplied by the length of the subsegment. This pro- 
cedure and the interpolation orders have been chosen in such a way that the 
viscous twisting and bending forces, which depend on the derivatives of order 
up to four of the positions, remain smooth upon subdivision. 



7. Interaction of the thread with other bodies 

Interaction with two types of rigid bodies are considered in our validation 
examples: a fluid container, such as a syringe, delivering fluid at a fixed volume 
rate, and a rigid impenetrable surface onto which the thread impinges. The 
interaction is taken care of using kinematical constraints. The vertices that have 
not yet been pushed outside the container, as well as those that have collided 
with the surface, are constrained. Their positions and velocities are prescribed 
by the motion of the body. In this section, we explain how the matrix i? and 
the vector B' , which represent these constraints in algorithm [l] are calculated. 
For the sake of simplicity, we assume that the motion of the external bodies is 
unaffected by the motion of the thread; two-way coupling with external bodies 
could be enforced as a kinematical constraint, and would typically be handled 
by a post-integration step using a manifold projection method |75j . 



7.1. Fluid container 

In the experiments, the viscous fluid is fed from a container at a constant 
volume rate Qc through a circular opening of diameter dc by a syringe controlled 
by a step motor, see figure [t^. Let = ^ be the area of the opening. The 
imposed volume rate sets an ejection velocity Uc = Qc/^c of the thread relative 
to the container. The velocity of the fluid exiting from the container is the sum 
of the prescribed velocity of the container, and this relative velocity Uc- 

At the exit from the container, the thread is clamped, i. e. the centerline 
position, tangent, and the rotation are all prescribed. This clamped boundary 
condition is implemented by blocking the rotation of the segment joining the 
first two vertices, and constraining the positions of these two vertices (this has 
the effect of blocking the unit tangent vector parallel to the first segment). These 
kinematical constraints are implemented as explained in section (6.2 1. Clamped 
boundary conditions therefore require that two vertices lie inside the container 
at all times. 



7.1.1. Simple container model 

A simple discrete model for the container is sketched in figure [7]d. The two 
topmost vertices of the thread are located inside the container and move with 
the velocity — C/cg^ relative to the container. Here is the unit vector, usually 
vertical, which is opposite to the direction of ejection. At every time step. 



49 



Figure 7: Discrete representation of the container delivering fluid through a circular opening of 
diameter dc at a prescribed vertical velocity Uc, shown in (a). Clamped boundary conditions 
are enforced by prescribing the position and velocities of the two topmost vertices, and blocking 
the rotation of the first segment joining them. Filled disks represent constrained vertices 
(top part of the thread attached to the container), open circles represents unconstrained 
vertices (hanging part of the thread), dashed curves are the timelines of vertices, (b) Simple 
implementation whereby the two first vertices move with prescribed vertical velocity Uc- (c) 
In a refined implementation, the first two vertices are kept at a fixed position with respect to 
the container, and mass is continuously added into the second segment. In both cases (b) and 
(c) a new vertex is periodically created from top (vertex labelled '0' appearing at time t + 2e), 
and there are always two vertices inside the container. 

whenever one of these two vertices goes past the opening of the container, it 
is freed (its position constraint is discarded) and a new constrained vertex is 
added on top of the thread, as happens in frame t + 2e in figure [7)d; the new 
segment on top is assigned a length 4, a volume {A^ic), a mass {pAcic), and 
surface tension 7. Here, p and 7 are the mass density and surface tension of the 
fluid in the container, and the length £c is a discretization parameter chosen by 
the user. 

This simple implementation has a drawback: it induces oscillations of the 
thread with a small amplitude and a large frequency. Indeed the effective fall 
height is determined by the position of the second topmost vertex. As a function 
of time, this height varies abruptly every time a new vertex is added. The 
amplitude of the discontinuity is set by the spacing between vertices in the 
containers. We found that this small-amplitude, large-frequency forcing usually 
allowed convergence of positions and velocities with finer and finer discretization 
lengths ic, but induces unbounded fluctuations of acceleration. 

7.1.2. Refined container model 

We found it necessary to suppress the oscillations to obtain convergent and 
fully reproducible results, and to capture the subtle patterns produced by the 
viscous sewing machine at large fall heights |76j . To do so, we used an improved 
description of the container which is sketched in figure [Tj:. The first two vertices 
do not follow fluid particles but instead stay at a fixed position with respect 
to the container. As earlier, their velocity is constrained to the value —UcC^ 
relative to the container at the start of each dynamic step, but their position is 
systematically reset after each dynamic step, so as to make the second vertex 
coincide with the opening of the container and the first vertex lie at distance 
ic above it. At each time step e, an incremental volume of fluid (eQc) and the 
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associated mass {peQc) provided by the container are assigned to the second 
segment. Whenever its volume exceeds the target value {Ac ic) fixed by the dis- 
cretization parameter ic, this segment is split, a new vertex is inserted between 
the second and third vertices; the fluid material in excess is then assigned to 
the new segment, as sketched in figure [7|;. In this improved implementation of 
the boundary conditions, the fall height varies smoothly as a function of time. 
We found that this was sufficient to make the acceleration converge smoothly in 
the limit of a fine discretization, and produce stitching patterns that are both 
reproducible and consistent with the sewing machine experiments. All examples 
shown in Sections [8] and |9] make use of this refined model. 

7.2. Collisions on a hard surface 

We now consider collisions with a rigid obstacle. When the thread hits the 
obstacle, it sticks to it and gets carried away by the obstacle. The obstacle may 
be at rest, as in the case of steady coiling of section[8j or in motion, as happens 
with the moving belt in the sewing machine experiment of section |9.2[ In the 
latter case, the motion of the obstacle is prescribed. Detection of collisions 
involves comparing the distance of vertices to the obstacle, to the radius a* of 



the adjacent segments defined in equation (81 1 



7.2.1. 'Capture and continue' mode 

In a straightforward implementation, called the 'capture and continue' mode, 
we detect collisions with the obstacle at the end of every dynamic step, and mark 
the colliding vertices as being captured: in all subsequent dynamic steps, their 
positions and velocities are constrained based on the motion of the obstacle 
and on their relative position to the obstacle when they were captured. Slip- 
free conditions are considered: the rotation of segments joining any pair of 
captured vertices is blocked as well. We found that this simple description of 
the collisions was another source of large, spurious fluctuations in acceleration. 
The oscillations can be interpreted by the fact that the vertical momentum 
resulting from a collision is not transferred to the thread until the following 
time step, when the position constraints take effect. Another drawback of the 
method, which is partly responsible for oscillations, is that the thread penetrates 
into the obstacle by a small but very irregular depth, roughly proportional to 
the time step duration. 

7.2.2. 'Time roll-back' mode 

The oscillations were removed by using a roll-back, which allows for a more 
accurate handling of collisions. A roll-back discards any dynamic time step 
ending up in unexpected collisions. The time step is recomputed with the motion 
of the colliding vertices constrained in such a way that they come exactly in 
contact with the obstacle at the end of the time step. Roll-back can be viewed 
as an iterative predictor-corrector or iterative constraint refinement method |77] : 
a new unexpected collision may take place during the second tentative time step, 
inducing a third attempt, etc. A list of expected collisions is kept and updated 
after unsuccessful time step. 
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Roll-back removes the two main limitations of the straightforward 'capture 
and continue' implementation: it suppress the delay in transferring momentum 
from the obstacle to the thread, and makes the thread land exactly at the sur- 
face of the obstacle, removing the unwanted rugosity produced by the 'capture 
and continue' mode. We found that roll-back does indeed suppress the spuri- 
ous oscillations very effectively, bringing about benefits similar to those of the 
improved nozzle implementation discussed in section 7.1.2 All the validation 
examples shown in the following sections make use of the time roll-back method. 

If many collisions occur during a single time step, the roll-back method 
does not work well and spurious fluctuations reemerge in the acceleration as a 
function a time. The reason is that we work in the context of a linearized implicit 
scheme, in which the equations of motion are always linearized using the solution 
at the start of the time step. The successive iterations of a given time step all 
start from the same configurations. As a result, they all make use of the same 
set of linearized equations. If too many collisions occur before the end of the 
time step, this linearization becomes a poor approximation of the actual motion. 
To work around this limitation, we combined roll-back with time adaptation: 
if a new collision is detected before the end of the tentative time step, the 
next tentative time step is shortened and scheduled to end at the collision time 
estimated from the previous iteration. This time adaptation reduces the time 
step dynamically, in such a way that there is at most one collision per time step. 
The drawback is the increased complexity in implementation, and the fact that 
the simulation time increases with the rate of collisions. 



8. Validation in a steady coiling geometry 

We proceed to validate our discrete model and verify our implementation, by 
checking convergence in the smooth limit. We consider the steady coiling motion 
of a viscous thread stretched by gravity and impinging on a surface at rest, as 
shown in figure [8}3. Our simulation results are compared to reference solutions 
kindly provided by N. Ribe, which are based on numerical continuation of the 
time-independent problem expressed in the co-rotating frame [9], and solved 
using the AUTO software [P^. 

8.1. Validation of bending, stretching, gravity, inertia and collisions 

The following set of parameters are used for validation and verification: the 
fluid's dynamical viscosity /i = 0.2 and mass density p = 5 10^^, the acceleration 
of gravity g = 9.81, the area Ac = 6.44 10"'^ of the circular outlet of the container 
and the fluid's volume rate Qc = 3.96 10^"^. The surface tension 7 is set to zero 
until we validate surface tension later in section [8.31 

Three dimensionless groups characterize the properties of the fluid and the 
container |TT : 
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where v /i/p is the kinematical viscosity; the diameter of the container's 
circular opening dc and the extrusion velocity Uc — Qc/Ac have been defined in 
section |7.1| The numerical values of the dimensionless groups are Hi = 7000, 
112 = 7 and Ha = for the set of parameters listed above. This corresponds to 
significant, but not extreme, stretching: the radius decreases by a factor of order 
2 during the course of the descent for the range of heights considered below. 

Based on the acceleration of gravity g and on the kinematical viscosity z/, 
one can define a natural length scale L* and a natural time scale T* by 

"' = (7) • ^' = t) ■ <™' 

Numerically, L* = 25.36 and T* — 1.61 using the above set of parameters. These 
scales are used to make the simulation results dimensionless when comparing to 
the reference solution. 

Two additional discretization parameters are needed in the simulation: the 
initial segment length ic, introduced in section [7.11 and the time step e. Unless 
otherwise specified, their values are set to ic = 0.025 and e = 0.02. Note that the 
average number of particles emitted per time step is (eUc/ic) — 0.49. A good 
trade-off between accuracy and efficiency requires that this number is neither 
very small nor very large. 

The reference solution of N. Ribe [5) takes the form of data for the coiling 
radius or the coiling frequency, as a function of the fall height. To produce 
simulation data which can be compared to this reference solution, a range of 
fall heights H(t) is swept in a single simulation run. To this end, the motion of 
the container is prescribed in a sequence of up to three phases in our simulations. 
At initial time t = 0, it is placed at height Hq. It is left still until time t — ti, 
when the steady coiling is established. Then, the container is moved upwards 
at a prescribed velocity Vc that is much smaller than the extrusion velocity Uc-, 
until time t = t2. In a last phase, from time t — t2to the end of the simulation, 
t = tf, the container is moved slowly downwards at the same velocity Vc- Unless 
otherwise specified, we set ti — 30, ^2 — 1530, tf — 3030 and Vc — 0.02, which 
is 30 times slower than the extrusion velocity. 

Validation results are shown in figure|8] Most of the time, the simulation lays 
down a thin curve (brown and blue) in the plane {H, R), where R = (x^ + y^)^/^ 
is the distance of the point of contact of the thread with the floor, to the axis 
passing through the nozzle. This indicates a steady coiling regime. When the 
container is moved up (brown curve) , the simulated radius R follows closely the 
reference solution, until the latter folds back onto itself. The portion of the 
reference curve immediately past the fold point is known to be unstable |llj . 
The simulation then goes to a transient regime, labeled 'C in figure |8^ and 
shown in figure [TT] It then settles to a different branch of solutions having a 
smaller radius. A series of such bifurcations is observed in the simulation, a 
behavior that has been observed in experiments too [11] . A similar sequence of 
transitions is observed when the container is moved down. As multiple coiling 
solutions are in competition, a hysteretic behavior is observed: transition from 
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Figure 8: Validation in a steady coiling geometry, with Hi = 7000, 112 = 7 and Ha = (no 
surface tension), (a) The coiling radius R = (x^ + y'^)^^'^ is recorded continuously as the fall 
height H is varied, by slowly moving the container upwards (brown curve) and then downwards 
(blue curve). It is compared against the reference solution of N. Ribe obtained by numerical 
continuation (black), after rescaling by the length scale L*. Note the abrupt changes in the 
coiling radius at places where the reference solution folds back onto itself, as expected. Filled 
regions correspond to a rapidly varying radius, either because the thread goes to a steady 
folding mode, as happens around H = .60, or because of a transient regime following a branch 
jump, (b) Two typical configurations of the thread, for different fall heights. Note that the 
fundamental mode of a hanging viscous string is excited in 'A', while its first harmonic is 
excited in 'B': the red curve has a node near z/L* = .05. 



the branch corresponding to the largest coiUng radius to the second largest 
occurs at a height 'C when the container moves up, which is larger than the 
transition height 'D' observed on the way back. We observe the occurrence of 
a folding mode in the interval .72 > H/L* > .47 when the container moves 
down, but will not comment further on it as little is known on the competition 
between the coiling and folding modes. The small gap between the blue and 
brown curves in the left part of figure [8^ can be attributed to the fact that 
the velocity of the container Vc is small but finite. Overall, the simulation is 
in good agreement with the reference curve, and reproduces the details of its 
meandering shape. This validates the various physical ingredients that affect 
the reference curve, namely viscous bending and stretching, gravity, inertia and 
contact with the floor. 

A detailed run such as the one shown in flgure [8] corresponding to a total 
simulation time — 3030 and a total number of time steps tf/e~ 150 10'^, runs 
in about 30 min on a 2.6 Ghz Intel Core 17 processor using 8 GB of memory. 
The maximum number of vertices, when the fall height is maximum, is 460, 
corresponding to approximately 1800 degrees of freedom. 

8.2. Analysis of convergence 

Convergence of the solution towards the reference solution of N. Ribe is 
shown in flgure [9] as a function of the discretization parameters. Convergence is 
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Figure 9: Analysis of convergence towards the reference solution in steady coiling geometry, 
showing the relative error on the coiling radius SR/R as a function of mesh size and time step. 
The same parameters are used as in figure |8] but the time step and the discretization length 
are varied in the intervals .02 < e < .2 and .025 < £c < -25, their ratio being kept constant. 
Fall height is fixed to H = 1.01. 




challenging in the presence of collisions on the floor, and we found it necessary 
to use the refined models for the floor and container described in sections 17.1.21 
and |7.2.2| In our simulation the coiling radius R is measured after a time that 
is long enough for the initial transient to disappear, and then averaged over 
several periods. The simulation is repeated for different values of the time step 
e and discretization length £c, their ratio being kept constant. Convergence of 
our numerical method is confirmed by the fact that the residual error goes to 
zero. The convergence appears to be linear with the discretization parameters. 

8.3. Validation of surface tension 

With the aim to validate our discrete model for surface tension, we repeat 
the validation shown in figure |8] using the same set of parameters, except for 
the surface tension coefficient, now set to 7 = 10~^. The corresponding dimen- 
sionless parameter is II3 — 10.310-"^. Surface tension has a marked effect on 
the coiling radius, as shown by comparison of the solid and dashed black curves 
in figure [lO) We obtain a good agreement between the simulation and the new 
reference curve. This validates our discrete surface tension model. 



9. Discussion 

9.1. Transient regimes 

Even though the steady coiling geometry provides a convenient set-up for 
validation and verification, our numerical method can solve the non-steady dy- 
namics of viscous threads. As an illustration, transient regimes following jumps 



from one branch of steady coiling solutions to the next are shown in figure 11 
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Figure 10: Validation of surface tension using steady coiling (Hi = 7000, 112 = 7, TI3 = 
10.3 10~^). The same parameters are used as in figure [s] except for the non-zero surface 



tension 7 



10" 



A good agreement is obtained with the reference curve that takes into 



account surface tension (solid black curve), 
surface tension is shown (dotted curve). 



For reference, the reference solution with zero 
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Figure 11: Trace laid out by the thread in the simulation, in the transient regimes following 
the jumps labelled 'C and 'D' in figure [s] In 'C, the fall height is slowly increased and the 
thread jumps to a solution having a smaller radius; in 'D', the fall height is slowly decreased, 
and the thread jumps to a solution having a larger radius. Same simulation parameters as in 
figure [s] 
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meanders alternating loops coiling 



Figure 12: Simulation of the so-called numerical sewing machine. The thread is poured from 
a constant height onto a moving belt. As the velocity of the belt is increased, the pattern laid 
out by the thread undergoes a series of bifurcations, similar to those that have been reported 
in the experiments: for this particular value of the fall height, translated coiling, alternated 
loops and meanders are successively obtained. Arrows indicate the velocity of the belt as the 
pattern was formed: the belt moves to the right and more ancient patterns, corresponding to 
a slower belt, are located on the right hand side. Simulation parameters are the same as in 
reference 1131 and are provided in main text. 



9.2. The viscous sewing machine 

The fluid-mechanical sewing machine is an extension of the steady coiling 
problem to the case of a moving substrate. To the best of our knowledge, our nu- 
merical method is the first one that can simulate this non-steady phenomenon. 
In fact the possibility to set up numerical simulations of the experiments de- 
scribed in references [12j [13] acted as an incentive for us to develop the present 
simulation method. 

A numerical investigation of the fluid-mechanical sewing machine is the topic 
of a separate detailed paper |76| . In figure 12 a simple demonstration run is 
presented. In this simulation, the velocity of the belt is steadily increased. 
The patterns must be read from right to left: this corresponds to the order in 
which they were produced, and so to an increasing belt velocities. On the right- 
hand side, for a low belt velocity, the translated coiling pattern is obtained. 
Increasing the belt velocity, two successive transitions are observed, leading to 
the formation of alternated loops first, and to meanders next. For even larger 
belt velocities, the oscillations disappear and the pattern becomes straight (not 
shown) as the hanging part of the thread takes on a catenary-like shape. The 
same sequence of patterns has been observed in the experiments, and is typical 
of small fall heights. A variety of patterns, some of which are quite complex, are 
obtained at larger fall heights both in the experiments and in the simulations, 
see reference (76i for details. 

The parameters used to produce the results of figure[T2]are the same as in the 
experiments of Morris [T3| . We use a convenient set of units in which the fluid's 
dynamical viscosity is /Lt = 1, its mass density is p ~ 1, and the acceleration 
of gravity is g = 1. With this natural set of units, the quantities L* and T* 
defined earlier read L* — 1 and T* = 1. The fall height is fixed to H = 0.865; 
this corresponds to a physical fall height of 3.7 cm in the experiments [13]. The 
area of the circular outlet at the bottom of the container is set to Ac = 0.0275, 
and the imposed volume flow rate is Qc — 2.2910^^, and there is no surface 
tension, 7 = 0. The corresponding values of the dimensionless groups read 
Hi = 608.5, 112 — 0.369, II3 — 0. The spatial discretization parameter is set to 
£c = 0.005, and the temporal one to e = 0.05. The stretching effect of gravity is 
severe, and the mesh is refined adaptively to maintain a good resolution near the 
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bottom of the thread while keeping a reasonable number of degrees of freedom, 



as explained in section 6.5 a given segment is refined up to six times in the 
simulation. The belt is held still until time ti — 919 to the initial transient 
relax. From this time on, the belt velocity is steadily accelerated until the end 
of the simulation, occurring at time tf = 1650, its final velocity being 0.02. 

9.3. Limitations and perspective 

For the simulation results to be reproducible, we found it critical to fine-tune 
the discretization parameters in a steady coiling geometry, before attempting to 
simulate non-steady behaviors. Checking the convergence of positions, typically 
by comparing the coiling radius to the reference solution, appeared not to be 
sufficient. We had to suppress any fast residual oscillation in the computed 
acceleration, by using the refined description of obstacles described in section [7] 
and by carefully choosing the discretization parameters. 

The present paper presents a linear implicit scheme: at each time step, 
the velocity is updated using a linear expression for the viscous forces, which 



is accurate in the vicinity of the previous configuration, see equations ( 126 



1271 and (135); the position is then updated separately using this velocity, see 
equation ( 136 ). Compared to an explicit scheme, this method is more difficult to 
implement but vastly superior in terms of stability. In a previous paper, we have 
explored a fully implicit scheme |41j , by retaining the non- linear dependence of 
the viscous forces on positions: in that case, each time step requires a non- 
linear root-finding. The benefit of the non-linear implicit approach is that it 
preserves conservation laws associated with the symmetries of the problem, such 
as the conservation of the angular momentum, when it applies. By contrast the 
method presented here displays the usual dissipation of angular momentum. 
This is a minor drawback as viscous thread are rarely free-standing in practical 
applications. In all the demonstration examples shown earlier, neither the linear 
nor the angular momentum of the thread are actually conserved, because of the 
contact forces with the floor. 

Simulations of viscous threads and elastic rods can be unified [41 : in that 
treatment all the code is shared by both rods and threads, except that in the 
viscous case the undeformed configuration is updated at each time step. This 
makes the combined approach very appealing. Nevertheless, a benefit of the 
specialized implementation presented here is that it saves the burden of imple- 
menting the Hessian matrix for naturally curved elastic rods. 

In future work, it would be interesting to extend the present numerical model 
to thin threads governed by more general constitutive laws, such as visco-elastic 
filaments |491 179| which can exhibit a complex and poorly understood behav- 
ior [S^. To this end, the discrete geometrical model exposed in the present 
paper can be reused and combined with different constitutive laws. It would 
also be interesting to couple our thin thread model with existing simulation 
methods for 3D flows with free boundaries [HI [82l [20l [72] , in order to capture 
the interaction of the thread with the slowly collapsing pile that forms where it 
merges with the bath. 
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9.4- Conclusion 

We have presented a numerical method for simulating the dynamics of thin 
viscous threads. In contrast with existing numerical methods, it captures the 
combined effects of stretching, bending and twisting forces, inertia and large 
rotations. It is not restricted to steady flows. The method has been derived by 
writing the smooth equations of motion for thin threads into a Lagrangian form, 
and using a careful spatial discretization. In particular, a discrete notion of twist 
has been used, which is based on the geometry of parallel transport. The internal 
stress representing the internal viscous stress has been derived from variational 
principles, using a Rayleigh potential. All the relevant physical quantities, such 
as strain rates and internal stress, have been identified in the discrete setting. 
Using a formal convergence argument, we have shown equivalence of the discrete 
equations with the classical smooth formulation. The method has been validated 
against reference solutions available for steady coiling. Demonstration examples 
in the non-steady case have been shown. 

We would like to thank Neil Ribe for getting us interested into the fascinating 
dynamics of thin threads, and for sharing his continuation data which enabled 
us to validate our code. 



Appendix A. Equivalence with the constitutive equations of Ribe 



In a classical paper, Ribe analyzed the helical coiling of viscous jets falling 
on a plane [9]. In the frame rotating with the jet, the shape of the centerline is 
stationary. The equations for the shape of the jet are expressed as a set of non- 
linear ordinary differential equations with boundary conditions at both ends. 
This non-linear boundary value problem was solved by numerical continuation 
techniques, using the AUTO j78j software. These solutions, which corresponds 
to steady configurations, have been used to validate our dynamical code. 

We show below that the constitutive laws derived by Ribe from the Stokes 
equations in 3D, are equivalent to our equations (65). However, the formalism 



used by Ribe is different from ours, and we then need to reword his analysis. 
Ribe introduces an Eulerian twist-curvature vector tt^ which satisfies an 



equation similar to our equation ( 28 ) for its Lagrangian variant tt, where the 



derivative is taken with respect to s instead of S, that is: 



ds 



(A.l) 



The components of tt^ are denoted 



Kl 
K3 



(A.2) 



where the square brackets indicate a decomposition in the material frame {di , d 
which is a 'moving' frame, i. e. this frame varies with the arc length parameter 
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As in equation ( 30 ) , the Eulerian kinematical twist and binormal curva- 
ture vector are obtained by a decomposition of the twist-curvature vector 
77^ into transverse and tangential components, 



K3, 







(A.3) 



In the frame rotating with the coils, the material rotation is written in a 
form similar to equation (27), 



where tf, is the time derivative evaluated in the rotating frame. Since the shape 
of the thread is steady in this frame, we have 



u 



ds 



Utt^ X t, 



where U = Mj^ • t is the axial velocity of the fluid in the rotating frame. The 
expression for the material angular velocity can be obtained by combining the 
above equations, 

U>3 



^R 



(A.4) 



The material angular velocity in the laboratory frame follows from the compo- 
sition of velocities, that is w = Wj^ + fle^ where fl is the frequency of coiling, 
as well as the relative angular velocity of the rotating frame with respect to the 
laboratory. 

As a result, the gradient of rotation reads 



did 

ds 



d_ 
ds 



' Uki' 








' Uki ' 








UK2 


Ll>3 




^3 




UJ3 



where the last term comes from the fact that the material frame is a moving 



frame, see equation (A.l). 

With the help of equations ( |A.2 ) and ( A.4[ ), this vector reads 



doj 
ds 



{U KlY + K2 {UJ3 -~ U K3) 

{U Ka)' - K3 {uj3 ~ U K3) 



When this vector is inserted into the constitutive law ( 65b I , we obtain for the 
bending moment m the same expressions as those derived by Ribe, who used a 
lubrication-type of approximation in the Stokes equations in 3D. 
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